- Email: [email protected]
Pseudo-rigid and rigid—flexible bodies
!i•iiii•ii••iiiii•i!••!i•ii•i!ii•i!iii!•ii!i•i•i••ii•ii••iii•••i! i~i i~i i~i i i ~i i~i!i i i i i !~ii i~!i~i~i i i~i !~i~i!~i~i~i i~i~i~i!i~i~i~
Pseudo-rigid and rigid-flexible bodies :,ii::~,iiiii''~i:,','i:,,iiii::i:,8i~ :::~:~' il::'~:::~:':iii!i!!'/' ii,i:::,i ,~,::i::ii'~:,'iii~, :i'~':i!''~::',i:~':i~',~"~i'~~,i~:',ii::~ii~ii ,i!'~,~,~ii~,:,~i,di'~i,:i'i~,:~:,:,i~~'i,:'~i,:,':ii,:i!i,Jii~i:~::il~~,'~,~i~,:i,'~!:,i',',i',i'i:,:~,i'::,':,:'i:,,i~:::iii i: i~'~!i:,:~:~':ii''i,i',!!',:,i:,, ,~i:'~::': ',':',i:,i'!,'~:~':,iii~:'~i'iii::i , :iii::::~i i~,i'':,:~i'i',:,::,~:~::i:~', ,i':!::i' ::::i'~':!::~::::~,i',i:,':~i:,~:':,~:::,:i:,~'~:, ,',~,ii','~::,'i,:,i~ii:,: i:i,::,::i:,:iii::~ii i ',:,i':,:,::'i~;,i:ii':,::,':i':',~::i:,:,::,~',i'i':i!'!i:,, !iiii', i:, :,~,il:,::::i:: : i ::i::ii: :,~ili::iii',i~'~~ii:,':,'~:, ::,iiii:,i',:,i!',:,i:,,~:,,~:~:,i:~, : ':'::,:'~','i:::i::ii i~:i::i:,ii::::~!i:ii:,ii',:,ii::,:,~,',',:i''~':,i::,!i::::i' '::,:~'::,iii:~:iil','i,:,!i'i':i',',i:il~:,',:!i,i':::,~i':i'',:if'i ::',i~:iil , :,i:,i:i::iiiiiiii:i:i:,,~:,'iliii:, '~i',:~,!:,i::':: :::ii::if:i~,::::i'ii:,i:'i::~::i:,!iilli'': i:::i'~ ,!:,::::,:~::i::ii' 'i: ::i'::i~iil
Many situations are encountered where treatment of the entire system as deformable bodies is neither necessary nor practical. For example, the frontal impact of a vehicle against a barrier requires a detailed modelling of the front part of the vehicle but the primary function of the engine and the rear part is to provide inertia, deformation being negligible for purposes of modelling the frontal impact. A second example, from geotechnical engineering, is the modelling of rock mass landslides or interaction between rocks on a conveyor belt where deformation of individual blocks is secondary. In this chapter we consider briefly the study of the first class of problems and in the next chapter the second type in much more detail. The above problem classes divide themselves into two further subclasses: one where it is necessary to include some simple mechanisms of deformation in each body (e.g. an individual rock piece) and the second in which the individual bodies have no deformation at all. The first class is calledpseudo-rigid body deformationl and the second rigid body behaviour. ~-Here we wish to illustrate how such behaviour can be described and combined in a finite element system. For the modelling of pseudo-rigid body analyses we follow closely the work of Cohen and Muncaster ~ and the numerical implementation proposed by Solberg and Papadopoulos. 3 The literature on rigid body analysis is extensive, and here we refer the reader to papers for additional details on methods and formulations beyond those covered here. 4-21
!iiiiiiiiiiiii~iiiii~iiiiiiiiiii~i~iiiiiiiiiiii~iiii~iii!!~i~iiiiiiiiiii~iiiiiii~iiii!ii!ii~iii
iiiiiiiiiiiiiiiiiiiii!~iii~iiiiiiiiiiiiiiiiiiiii!i!iiiiiiiiiiiiiiii!iiiiii!!iiiiiiiiiiiiiii i i iiiiii!i!iiiiiii i i iiiiiiiiiiiiiiiiii iiii iiiiiiiii!i iiii i!iiiiiiiii!ii!ii i !iiiiiii!iiiiii iiiiiiiiiiiiiiiiiiiiiiii i !ii
~T~,,~.~ v.,~L~',~~,,,,seuao-rlgm~motlOnS~'~,~
','~',~,~',~',~,~ ~~
~',~
;~ ~', ~ ~'~',~ ~ ~',~,~',~,~,~,~ ~ ~'~', ,~',
In this section we consider the analysis of systems which are composed of many small bodies, each of which is assumed to undergo large displacements and a uniform deformation.* The individual bodies which we consider are of the types shown in Fig. 8.1. In particular, a faceted shape can be constructed directly from a finite element discretization in which the elements are designated as all belonging to a single solid object or the individual bodies can be described by simple geometric forms such as discs or ellipsoids. * Higher-order body.
approximations
can
be
included
using
polynomial
approximation
for
the
deformation
of each
Pseudo-rigid motions 229
(a)
(b)
Fig. 8.1 Shapesfor pseudo-rigid and rigid body analysis: (a) ellipsoid; (b) faceted body. A homogeneous motion of a body may be written as
r
t) = ri(t) -+- Fii(t) [XI - RI]
(8.1)
in which XI is position, t is time, RI is some reference point in the undeformed body, ri is the position of the same point in the deformed body, and Fil is a constant deformation gradient. We note immediately that at time zero the deformation gradient is the identity tensor (matrix) and Eq. (8.1) becomes
r
O) - -
ri(O) -+- tSiI [ X I -- RI] = ri(O) + t ~ i l X l -- tSilRl ~ t S i l X l
(8.2)
where ri(O) "-- t~iiRi by definition. The behaviour of solids which obey the above description is sometimes referred to as analysis ofpseudo-rigid bodies. 1 A treatment by finite elements has been considered by Solberg and Papadopoulos, 3 and an alternative expression for motions restricted to incrementally linear behaviour has been developed by Shi, and the method is commonly called discontinuous deformation analysis (DDA). 22 The DDA form, while widely used in the geotechnical community, is usually combined with a simple linear elastic constitutive model and linear strain-displacement forms which can lead to large errors when finite rotations are encountered. Once the deformation gradient is computed, the procedures for analysis follow the methods described in Chapter 5. It is, of course, necessary to include the inertial term for each body in the analysis. No difficulties are encountered once a shape of each body is described and a constitutive model is introduced. For elastic behaviour it is not necessary to use a complicated model, and here use of the St Venant-Kirchhoff relation is adequate - indeed, if large deformations occur within an individual body the approximation of homogeneous deformation generally is not adequate to describe the solution. The primary difficulty for this class of problems is modelling the large number of interactions between bodies by contact phenomena and here the reader is referred to Chapters 7 and 9 and references on the subject for additional information on contact and other details. 22'23
230
Pseudo-rigidand rigid-flexible bodies
The pseudo-rigid body form can be directly extended to rigid bodies by using the polar decomposition on the deformation tensor. The polar decomposition of the deformation gradient may be given as 24-26 Fil -- A i j U j I
where
AilAij
"- (~lJ
and
AilAjl
= ~Sij
(8.3)
Here AiI is a rigid rotation* and UIj is a stretch tensor (which has eigenvalues Am as defined in Chapter 5). In the case of rigid motions the stretches are all unity and UIj simply becomes an identity. Thus, a rigid body motion may be specified as (gi(Xl, t) = ri(t) + A i l ( t ) [XI - RI]
(8.4)
q~(X, t) = r(t) + A ( t ) [ X - R]
(8.5)
or, in matrix form, as Alternatively, we can express the rigid motion using Eq. (8.1) and impose constraints to make the stretches unity. For example, in two dimensions we can represent the motion in terms of the displacements of the vertices of a triangle and apply constraints that the lengths of the triangle sides are unchanged during deformation. The constraints may be added as Lagrange multipliers or other constraint methods and the analysis may proceed directly from a standard finite element representation of the triangle. Such an approach has been used in reference 27 with a penalty method used to impose the constraints. Here we do not pursue this approach further and instead consider direct use of rigid body motions to construct the formulation. For subsequent use we note the form of the variation of a rigid motion and its incremental part. These may be expressed as A
3~b = 3r + 30A [X - R] dO - dr + dOA [X - R]
Using Eq. (8.5) these may be simplified to 3q~ - 3 r - ~30 dc~ - d r - ~dO
where
y = x- r (8.6)
where dO and ~0 are incremental and variational rotation vectors, respectively. In a similar manner we obtain the velocity for the rigid motion as -/"-
~w
(8.7)
in which f is translational velocity and w angular velocity, both at the centre of mass. The angular velocity is obtained by solving h = a,A * Often literature denotes this rotation as Rii; however, here we use to avoid confusion use Ail to denote rotation.
(8.8) RI
as a position of a point in the body and
Rigid motions 231 or
A-
A~
(8.9)
where f~ is the reference configuration angular velocity. 8 This is clearer by writing the equations in indicial form given by ~kil -- w i j A j l
(8.10)
-- A i j ~ ' 2 j l
where the velocity matrices are defined in terms of vector components and give the skew symmetric form Wij =
W3 --W 2
0
--W1
W1
0
(8.11)
and similarly for f2lj. The above form allows for the use of either the material angular velocity or the spatial one. Transformation between the two is easily performed since the rigid rotation must satisfy the orthogonality conditions ATA = A A T = I
(8.12)
at all times. Using Eqs (8.8) and (8.9) we obtain -- A ~ A T
(8.13)
= ATcDA
(8.14)
or by transforming in the opposite way
8.3.1 Equations of motion for a rigid body If we consider a single rigid body subjected to concentrated loads fa applied at points whose current position is Xa and locate the reference position for R at the centre of mass, the equations of equilibrium are given by conservation of linear momentum = ~
fa = f;
p = mr
(8.15)
a
where p defines a linear momentum, f is a resultant force and total mass of the body is computed from P
m - ]
J.
P0 dV
(8.16)
and conservation of angular momentum ~" - - Z
(Xa - -
r) • fa -- m;
Jr = I[w
(8.17)
a
where rr is the angular momentum of the rigid body, m is a resultant couple and ]I is the spatial inertia tensor.
232
Pseudo-rigidand rigid-flexible bodies The spatial inertia tensor (matrix) ITis computed from ]I = A J A T
(8.18)
where J is the inertia tensor (matrix) computed from an integral on the reference configuration and is given by ,~ - J~ P0 [(YTY) I -- yyT] dV
where
Y - X - R
(8.19)
Thus, description of an individual rigid body requires locating the centre of mass R and computing the total mass m and inertia matrix Jl. It is then necessary to integrate the equilibrium equations to define the position r and the orientation of the body A.
8.3.2 Construction from a finite element model If we model a body by finite elements, as described throughout this volume, we can define individual bodies or parts of bodies as being rigid. For each such body (or part of a body) it is then necessary to define the total mass, inertia matrix, and location of the centre of mass. This may be accomplished by computing the integrals given by Eqs (8.16) and (8.19) together with the relation to determine the centre of mass given by mR -- f n poX d V
(8.20)
In these expressions it is necessary only to define each point in the volume of an element by its reference position interpolation X. For solid (e.g. brick or tetrahedral) elements such interpolation is given by Eq. (5.50a) which in matrix form becomes (omitting the summation symbol) X = N.X~
(8.21)
This interpolation may be used to determine the volume element necessary to carry out all the integrals numerically by quadrature. 28 The total mass may now be computed as
m=~-'~(f~ podV) e
(8.22)
e
where ~"2 e is the reference volume of each element e. Use of Eq. (8.21) in Eq. (8.20) to determine the centre of mass now gives
R - --I _ m
e
PoN,~ dV) X~
(8.23)
e
and finally the reference inertia tensor (matrix) as J = E
M~;~ [(YSY~)I- Y~Y~] ; e
Y~ = X~ - R
(8.24)
Rigid motions 233 where
Mea3 = f~ PoN"N~dV
(8.25)
e
The above definition of Y, tacitly assumes that y ~ N, - 1. If other interpolations are used to define the shape functions (e.g. hierarchical shape functions) it is necessary to modify the above procedure to determine the mass and inertia matrix.
8.3.3 Transient solutions The integration of the translational rigid term r may be performed using any of the methods described in reference 28 or indeed by other methods described in the literature. The integration of the rotational part can also be performed by many schemes; however, it is important that updates of the rotation produce discrete time values for rigid rotations which retain an orthonormal character, that is, the An must satisfy the orthogonality condition given by Eq. (8.12). One procedure to obtain this is to assume that the angular velocity within a time increment is constant, being measured as
w(t) ~
(,~t3n+c~ =
•At
(8.26)
in which At is the time increment between tn and tn+l, 0 is the increment of rotation during the time step, and 0 < c~ < 1. The approximation (at,~nq_c~ --
(8.27)
(1 -- OL)(,~3 n -'1- OL(aL3n+ 1
is used to define intermediate values in terms of those at tn and tn+l. Equation (8.8) now becomes a constant coefficient ordinary differential equation which may be integrated exactly, yielding the solution A(t) -- exp[0(t - tn)/At]An
tn < t < tn-t-1
(8.28)
In particular at tn+~ we obtain An+c~ - exp[ce0]An This may also be performed using the material angular velocity f~.8 Many algorithms exist to construct the exponential of a matrix, and the closed-form expression given by the classical formula of Euler and Rodrigues (e.g. see Wittaker 29) is quite popular. This is given by
exp[0] - I +
1010 + _1 sin 2 101/202 101 2 [1012/2]
sin
where
101 =
[0T0 ] 1/2
(8.29)
This update may also be given in terms of quaternions and has been used for integration of both rigid body motions as well as for the integration of the rotations appearing in three-dimensional beam formulations (see Chapter 17). 8'30'31Another alternative to the
234
Pseudo-rigidand rigid-flexible bodies direct use of the exponential update is to use the Cayley transform to perform updates for A which remain orthonormal. Once the form for the update of the rigid rotation is defined any of the integration procedures defined in reference 28 may be used to advance the incremental rotation by noting that 0 or 19 (the material counterpart) are in fact the change from time t~ to tn+~. The reader also is referred to reference 8 for additional algorithms directly based on the GN 11 and GN22 methods. 28 Here forms for conservation of linear and angular momentum are of particular importance.
In some analyses the rigid body is directly attached to flexible body parts of the problem [Fig. 8.2(a)]. Consider a rigid body that occupies the part of the domain denoted a s ~2 r and is 'bonded' to a flexible body with domain f2f. In such a case the formulation to 'bond' the surface may be performed in a concise manner using Lagrange multiplier constraints. We shall find that these multiplier constraints can be easily eliminated from the analysis by a local solution process, as opposed to the need to carry them to the global solution arrays as was the case in their use in contact problems (see Sec. 5.3).
8.4.1 Lagrange multiplier constraints A simple two-dimensional rigid-flexible body problem is shown in Fig. 8.2(a) in which the interface will involve only three-nodal points. In Fig. 8.2(b) we show an exploded view between the rigid body and one of the elements which lies along the rigid-flexible interface. Here we need to enforce that the position of the two interface nodes for the element will have the same deformed position as the corresponding point on the rigid body. Such a constraint can easily be written using Eq. (8.4) as Ca = r(t) + A(t) [X~ - R] - x~(t) = 0
(8.30)
in which the subscript c~ denotes a node number. We can now modify a functional to include the constraint using a classical Lagrange multiplier approach in which we add the term 1-Irf = )~C~ = ,,~ [x~(t) - r(t) - A(t) [X~ - R]]
(8.31)
0
(a)
(b)
Fig. 8.2 Lagrange multiplier constraint between flexible and rigid bodies: (a)rigid-flexible body; (b) Lagrange multipliers.
Connecting a rigid body to a flexible body 235 Taking the variation we obtain
aFlrf = aA~ [x~ - r - A [X~ - R]] + A~ [ax~ - 6r - a0A [X~ - R]]
(8.32)
From this we immediately obtain the constraint equation and a modification to the equilibrium equations for each flexible node and the rigid body. Accordingly, the modified variational principle may now be written for a typical node a on the interface of the rigid body as ~I'I "4- ~I'Irf :
[~Xtz
6x~
~r
dO ~X,]
M.~v~ + M ~ v ~ + P. - f~, M ~ v ~ + M~av/~ + P~ - f~ + X~ (8.33)
#-f-X~
=0
e-m-yTA x~ - r -
A [ X ~ - R]
where y~ = x~ - r are the nodal values of y, fi are any other rigid body nodes connected to node a and #, v are flexible nodes connected to node a. Since the parameters x~ enter the equations in a linear manner we can use the constraint equation to eliminate their appearance in the equations. Accordingly, from the variation of the constraint equation we may write
6x~ = [I, y~]
(8.34)
60
which permits the remaining equations in Eq. (8.33) to be rewritten as
Mv~v~ + Mv~'b~ + P~ - f~ ~1"1 + ~I-Irf -- [~Xtt
&
~0]
p -- f + M~jr
+ M~av~ + P~ - f~
= 0
"/r - m - ~,T(M,~,% + M,~zi'a + P~ - f~)
(8.35) For use in a Newton solution scheme it is necessary to linearize Eq. (8.35). This is easily achieved
d(~I'I)+d(~I-Irf)= [ax~
ar I
30]
(K~)T
(Ktzfl)T
~X~
0
(Kav)T (Ka/~)T K.~ --yT(Kav)T --~T(Ka/~)T 0
qx~ KTo
dr
dO
(8.36)
236 Pseudo-rigidand rigid-flexible bodies Once again this form may be reduced using the equivalent of Eq. (8.34) for an incremental dx~ to obtain
d(~Fl) + d(~l-Irf )
[~x,
8r ~0]
r
(K~)T --(Kvr162 (KIzv)T - (K,~)T~'~ / (K~v)T K.~ + (K~3)T ,,T K [K~ + L--YT(K~.)T -Y~( ~;~)T
x x
l |
(8.37)
dr
dO Combining all the steps we obtain the set of equations for each rigid body as I
-- (Ktzfl)T~r/3
(Kvv)T
(Kv~)T
(K~)T
[K~ + (K~)T]
_~,T (Kav)T
_~,T (Ka~)T
]fdxv 1 ~d;
--(Ka/~)T~r/3 / [K~ + :9~(K~)T:9~IJ
(8.38)
=
m - ~ + ~'T'I'~ OL
in which ~ and ~ . are the residuals from the finite element calculation at node a and #, respectively. We recall from Chapter 5 that each is given by a form (8.39) which is now not zero since total balance of momentum includes the addition of the
~a. The above steps to compute the residual and the tangent can be performed in each element separately by noting that
,,ks = ~
Xe
(8.40)
e
where Xe denotes the contribution from element e. Thus, the steps to constrain a flexible body to a rigid body are once again a standard finite element assembly process and may easily be incorporated into a solution system. The above discussion has considered the connection between a rigid body and a body which is modelled using solid finite elements (e.g. quadrilateral and hexahedral elements in two and three dimensions, respectively). It is also possible directly to connect beam elements which have nodal parameters of translation and rotation. This is easily performed if the rotation parameters of the beam are also defined in terms of the rigid rotation A. In this case one merely transforms the rotation to be defined relative to the reference description of the rigid body rotation and assembles the result directly into the rotation terms of the rigid body. If one uses a rotation for both the beam and the rigid body which is defined in terms of the global Cartesian reference configuration no transformation is required. Shells can be similarly treated; however, it is best then to define the shell directly in terms of three rotation parameters instead of only two at points where connection is to be performed. 32'33
Multibody coupling by joints
Often it is desirable to have two (or more) rigid bodies connected in some specified manner. For example, in Fig. 8.3 we show a disc connected to an arm. Both are treated as rigid bodies but it is desired to have the disc connected to the arm in such a way that it can rotate freely about the axis normal to the page. This type of motion is characteristic of many rotating machine connections and it as well as many other types of connections are encountered in the study of rigid body motions. 4'34 This type of interconnection is commonly referred to as a joint. In quite general terms joints may be constructed by a combination of two types of simple constraints: translational constraints and rotational constraints.
8.5.1 Translation constraints The simplest type of joint is a spherical connection in which one body may freely rotate around the other but relative translation is prevented. Such a situation is shown in Fig. 8.3 where it is evident the spinning disc must stay attached to the rigid arm at its axle. Thus it may not translate relative to the arm in any direction (additional constraints are necessary to ensure it rotates only about the one axis - these are discussed in Sec. 8.5.2). If a full translation constraint is imposed a simple relation may be introduced as C j - - x (a) - x (b) - - 0 (8.41) where a and b denote two rigid bodies. Thus, addition of the Lagrange multiplier constraint (8.42) l - [ j - - ,,~jT[x(a) -- X(b)] imposes the spherical joint condition. It is necessary only to define the location for the spherical joint in the reference configuration. Denoting this as Xj (which is common
Fig. 8.3 Spinning disc constrained by a joint to a rigid arm.
237
238
Pseudo-rigidand rigid-flexible bodies to the two bodies) and introducing the rigid motion yields a constraint in terms of the rigid body positions as l-lj -- )kT[r (a) + A ( a ) ( x j
J
- R (a)) - r (b) - A(b)(xj -- R(b))]
(8.43)
The variation and subsequent linearization of this relation yields the contribution to the residual and tangent matrix for each body, respectively. This is easily performed using relations given above and is left as an exercise for the reader. If the translation constraint is restricted to be in one direction with respect to, say, body a it is necessary to track this direction and write the constraint accordingly. To accomplish this the specific direction of the body a in the reference configuration is required. This may be computed by defining two points in space X l and X2 from which a unit vector V is defined by V=
X2-Xl
IXz - Xll
(8.44)
The direction of this vector in the current configuration, v, may be obtained using the rigid rotation for body a v = A(a)v (8.45) A constraint can now be introduced into the variational problem as
I-lj =/kj{VT(A(a))T[r (a) + A(a)(xj - R (a)) - r O) - A~
- R(b))]}
(8.46)
where, owing to the fact there is only a single constraint direction, the Lagrange multiplier is a scalar/~j and, again, Xj denotes the reference position where the constraint is imposed. The above constraints may also be imposed by using a penalty function. The most direct form is to perturb each Lagrange multiplier form by a penalty tenn. Accordingly, for each constraint we write the variational problem as 1
2kj AJ
~
1-Ij -- )~jCj
2
(8.47)
where it is immediately obvious that the limit kj ---> oo yields exact satisfaction of the constraint. Use of a large kj and variation with respect to )~j give
E 1 ]
~%j C j - ~ ) ~ j
=0
(8.48)
and may easily be solved for the Lagrange multiplier as
)~j = kjCj
(8.49)
which when substituted back into Eq. (8.47) gives the classical form
nj -
kj
[cj
]2
(8.50)
The reader will recognize that Eq. (8.47) is a mixed problem, whereas Eq. (8.50) is irreducible. An augmented Lagrangian form is also possible following the procedures described in reference 28 and used in Chapter 7 for contact problems.
Multibody coupling by joints 239
8.5.2 Rotation constraints A second kind of constraint that needs to be considered relates to rotations. We have already observed in Fig. 8.3 that the disc is free to rotate around only one axis. Accordingly, constraints must be imposed which limit this type of motion. This may be accomplished by constructing an orthogonal set of unit vectors V I in the reference configuration and tracking the orientation of the deformed set of axes for each body as V~ c) --
~ilA(C)Vl
for c = a, b v~Vj - 5ij
(8.51)
A rotational constraint which imposes that axis i of body a remains perpendicular to axis j of body b may then be written as (v!a))Tv t
(b) - j
vT(A(a))TA(b)vj
-- 0
(8.52)
Example 8.1: Revolute joint As an example, consider the situation shown for the disc in Fig. 8.3 and define the axis of rotation in the reference configuration by the Cartesian unit vectors EI (i.e. V I = EI). If we let the disc be body a and the arm body b the set of constraints can be written as (where v3 is axis of rotation) X (a) __ X (c)
Cj =
~,*l[w(a)~Tw(b)J "3 = 0
(8.53)
l'w(a)~Tw(b) k*2 I *3
and included in a formulation using a Lagrange multiplier form
rlj = )~jT Cj
(8.54)
The modifications to the finite element equations are obtained by appending the variation and linearization of Eq. (8.54) to the usual equilibrium equations. Here five Lagrange multipliers are involved to impose the three translational constraints (spherical joint) and the angle constraints for the rotating disc. The set of constraints is known as a revolute joint, e
8.5.3 Library of joints Translational and rotational constraints may be combined in many forms to develop different types of constraints between rigid bodies. For the development it is necessary to have only the three types of constraints described above. Namely, the spherical joint, a single translational constraint, and a single rotational constraint. Once these are available it is possible to combine them to form classical constraint joints and here the reader is referred to the literature for the many kinds commonly encountered. 2,4,7,35 The only situation that requires special mention is the case when a series of rigid bodies is connected together to form a closed loop. In this case the method given above can lead to situations in which some of the joints are redundant. Using Lagrange multipliers this implies the resulting tangent matrix will be singular and, thus, one
240
Pseudo-rigid and rigid-flexible bodies cannot obtain solutions. Here a penalty method provides a viable method to circumvent this problem. The penalty method introduces elastic deformation in the joints and in this way removes the singular problem. If necessary an augmented Lagrangian method can be used to keep the deformation in the joint within required small tolerances. An alternative to this is to extract the closed loop rigid equations from the problem and use singular valued decomposition 36 to identify the redundant equations. These may then be removed by constructing a pseudo-inverse for the tangent matrix of the closed loop. This method has been used successfully by Chen to solve single loop problems. 35
'iii!iiiiiii iiiiiiiiiiii'ii,i'~i~ ,i~:~:',~::~::~iiiiiii'i,~:~i : :::~~i~,i,i:~,':ii':~iiii'i',i{iiii::~::::~'~'i~ii'~,iiiii::::~i~i'i~il~,ii,i~','i~,'i',i~',il!ili'~,i',iliiiiiii! i i:,i:,iiiiiiii~,~,',i~i{i'i,{!'!,!~,i ~!',iiiiiiiiii 'i,i ~""i~',~{'~,'!i~~:~ !i'i:~i',~'i~i'i i'~,i'~i,iiii'i!i'i,'i ':,,i'i,~i~,'il~,i!iiiiliiiiiiii' i i~i~,i',i!'i,:'iii'',',:iii!iii~'~'':iiili',iii',ii!'iii'~i',i',!i~,:,il,i',:{:i''i,~,i,i:ii!:i'',:,iii",','i, iii'!~,:',:ii',ii'"i~,,'~ili:,l'i,'ii'i:{~iiiiii'','i:,,',i!!,i~,i:i~:,'i:':,i:ii:i:,:i:,i::iif'~,i,'!',!i!'ii~'i,i'::,'i:~',,i'::',!,!iii!ii{i:i:,i!iii:,!:,i!':!:i{::i','i:i"~,ii'iiii!':i'~,'~,:!,i~:,i':i',i:i',,i','!'~,,i'i:!i~,,',',ii:: i"i'':,~,i i'!i,'','!~ii',i:,i' : ii',i'::,{:iiii'i':,,',i!:,:,i:,ii:!ii:i:' '~!~,i::ii':i':~ !i~i:'i :i,i','':iii',i,,'~~!: ':ii': i'
8.6.1 Rotating disc As a first example we consider a problem for the rotating disc on a rigid arm which is attached to a deformable base as shown in Fig. 8.4. The finite element model is constructed from four-node displacement elements in which a St Venant-Kirchhoff material model is used for the elastic part. The elastic properties in the model are E = 10 000 and z~ = 0.25, with a uniform mass density P0 = 5 throughout. The disc and arm are made rigid by using the procedures described in this chapter. The disc is attached to the arm by means of a revolute joint with the constraints imposed using the Lagrange multiplier method. The rigid arm is constrained to the elastic support by using the local Lagrange multiplier method described in Sec. 8.4. The problem is excited by a constant vertical load applied at the revolute joint and a torque applied to spin the disc. Each load is applied for the first 10 units of time. The mesh and configuration are shown in Fig. 8.4(a). Deformed positions of the model are shown at 2.5 unit intervals of time in Fig. 8.4(b)-(h). A marker element shows the position of the rotating disc. The displacements at the revolute joint and the radial exterior point at the marker element location are shown in Fig. 8.5.
8.6.2 Beam with attached mass ....................................................................................
--: .....................................................................................................................................................................................
As a second example we consider an elastic cantilever beam with an attached end mass of rectangular shape. The beam is excited by a horizontal load applied at the top as a triangular pulse for two units of time. The rigid mass is attached to the top of the beam by using the Lagrange multiplier method described in Sec. 8.4 and here it is necessary to constrain both the translation and the rotation parameters of the beam. The beam is three dimensional and has an elastic modulus of E = 100 000 and a moment of inertia in both directions of I1~ = 122 = 12. The beam mass density is low, with a value of /90 = 0.02. The tip mass is a cube with side lengths 4 units and mass density/90 = 1. The shape of the beam at several instants of time is shown in Fig. 8.6 and it is clear that large translation and rotation is occurring and also that the rigid block is correctly following a constrained rigid body motion.
Numerical examples 241
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Fig. 8.4 Rigid-flexible model for spinning disc: (a) problem definition, solutions at time; (b) t = 2.5 units; (c) t - 5.0 units; (d) t - 7.5 units; (e) t -- 10.0 units; (f) t - 12.5 units; (g) t - 15.0 units; (h) t - 17.5 units.
242
Pseudo-rigid and rigid-flexible bodies
15, -,r
0.5
E o
0
r 03
-0.5
A
A
l!
h
1 I I---Horizontall 10~"'Vertical I
A
".......",..-'I .- ::, ,,:I,- 1 0
4_~,,-
L
---Vertical J _ . -1 10 20 0 Time (a)
-15 -20 30
40
/""%
-.
'
0
10
'
20 Time (b)
30
40
Fig. 8.5 Displacements for rigid-flexible model for spinning disc. Displacement at: (a) revolute; (b) disc rim.
(a)
(c)
(b)
(e)
(f)
(g)
(d)
(h)
(i)
Fig. 8.6 Cantilever with tip mass: (a) t -- 2 units; (b) t -- 4 units; (c) t -- 6 units; (d) t - 10 units; (e) t -- 12 units; (f) t - 14 units; (g) t - 16 units; (h) t - 18 units; (i) t - 20 units.
1. H. Cohen and R.G. Muncaster. The Theory of Pseudo-rigid Bodies. Springer, New York, 1988. 2. A.A. Shabana. Dynamics of Multibody Systems. John Wiley & Sons, New York, 1989. 3. J.M. Solberg and P. Papadopoulos. A simple finite element-based framework for the analysis of elastic pseudo-rigid bodies. International Journal for Numerical Methods in Engineering, 45:1297-1314, 1999. 4. D.J. Benson and J.O. Hallquist. A simple rigid body algorithm for structural dynamics programs. International Journal for Numerical Methods in Engineering, 22:723-749, 1986.
References 243 5. R.A. Wehage and E.J. Haug. Generalized coordinate partitioning for dimension reduction in analysis of constrained dynamic systems. Journal of Mechanical Design, 104:247-255, 1982. 6. A. Cardona and M. Geradin. Beam finite element nonlinear theory with finite rotations. International Journal for Numerical Methods in Engineering, 26:2403-2438, 1988. 7. A. Cardona, M. Geradin and D.B. Doan. Rigid and flexible joint modelling in multibody dynamics using finite elements. Computer Methods in Applied Mechanics and Engineering, 89:395418,91. 8. J.C. Simo and K. Wong. Unconditionally stable algorithms for rigid body dynamics that exactly conserve energy and momentum. International Journal for Numerical Methods in Engineering, 31:19-52, 1991. [Addendum: 33:1321-1323, (1992).] 9. H.T. Clark and D.S. Kang. Application of penalty constraints for multibody dynamics of large space structures. Advances in the Astronautical Sciences, 79:511-530, 1992. 10. G.M. Hulbert. Explicit momentum conserving algorithms for rigid body dynamics. Computers and Structures, 44:1291-1303, 1992. 11. M. Geradin, D.B. Doan and I. Klapka. MECANO: a finite element software for flexible multibody analysis. Vehicle System Dynamics, 22:87-90, 1993. Supplement issue. 12. S.N. Atluri and A. Cazzani. Rotations in computational solid mechanics. Archives of Computational Methods in Engineering, 2:49-138, 1995. 13. O.A. Bauchau, G. Damilano and N.J. Theron. Numerical integration of nonlinear elastic multibody systems. International Journal for Numerical Methods in Engineering, 38:2727-2751, 1995. 14. J.A.C. Ambr6sio. Dynamics of structures undergoing gross motion and nonlinear deformations: a multibody approach. Computers and Structures, 59:1001-1012, 1996. 15. R.L. Huston. Multibody dynamics since 1990. Applied Mechanics Reviews, 49:$35-$40, 1996. 16. O.A. Bauchau and N.J. Theron. Energy decaying scheme for non-linear beam models. Computer Methods in Applied Mechanics and Engineering, 134:37-56, 1996. 17. O.A. Bauchau and N.J. Theron. Energy decaying scheme for non-linear elastic multi-body systems. Computers and Structures, 59:317-331, 1996. 18. C. Bottasso and M. Borri. Energy preserving/decaying schemes for nonlinear beam dynamics using the helicoidal approximation. Computer Methods in Applied Mechanics and Engineering, 143:393-415, 1997. 19. O.A. Bauchau. Computational schemes for flexible, nonlinear multi-body systems. Multibody System Dynamics, 2:169-222, 1998. 20. O.A. Bauchau and C.L. Bottasso. On the design of energy preserving and decaying schemes for flexible nonlinear multi-body systems. Computer Methods in Applied Mechanics and Engineering, 169:61-79, 1999. 21. O.A. Bauchau and T. Joo. Computational schemes for non-linear elasto-dynamics. International Journal for Numerical Methods in Engineering, 45:693-719, 1999. 22. G.-H. Shi. Block System Modelling by Discontinuous Deformation Analysis. Computational Mechanics Publications, Southampton, 1993. 23. E.G. Petocz. Formulation and analysis of stable time-stepping algorithms for contact problems. PhD thesis, Department of Mechanical Engineering, Stanford University, Stanford, California, 1998. 24. M.E. Gurtin. An Introduction to Continuum Mechanics. Academic Press, New York, 1981. 25. L.E. Malvern. Introduction to the Mechanics of a Continuous Medium. Prentice-Hall, Englewood Cliffs, NJ, 1969. 26. J. Bonet and R.D. Wood. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press, Cambridge, 1997. ISBN 0-521-57272-X. 27. J.C. Garcia Orden and J.M. Goicolea. Dynamic analysis of rigid and deformable multibody systems with penalty methods and energy-momentum schemes. Computer Methods in Applied Mechanics and Engineering, 188:789-804, 2000.
244 Pseudo-rigidand rigid-flexible bodies 28. O.C. Zienkiewicz, R.L. Taylor and J.Z. Zhu. The Finite Element Method: Its Basis and Fundamentals. Butterworth-Heinemann, Oxford, 6th edition, 2005. 29. E.T. Wittaker. A Treatise on Analytical Dynamics. Dover Publications, New York, 1944. 30. J.H. Argyris and D.W. Scharpf. Finite elements in time and space. Nuclear Engineering and Design, 10:456--469, 1969. 31. A. Ibrahimbegovic and M. A1 Mikdad. Finite rotations in dynamics of beams and implicit timestepping schemes. International Journal for Numerical Methods in Engineering, 41:781-814, 1998. 32. J.C. Simo. On a stress resultant geometrically exact shell model. Part VII: Shell intersections with 5/6 DOF finite element formulations. Computer Methods in Applied Mechanics and Engineering, 108:319-339, 1993. 33. P. Betsch, F. Gruttmann and E. Stein. A 4-node finite shell element for the implementation of general hyperelastic 3d-elasticity at finite strains. Computer Methods in Applied Mechanics and Engineering, 130:57-79, 1996. 34. H. Goldstein. Classical Mechanics. Addison-Wesley, Reading, 2nd edition, 1980. 35. A.J. Chen. Energy-momentum conserving methods for three dimensional dynamic nonlinear multibody systems. PhD thesis, Department of Mechanical Engineering, Stanford University, Stanford, California, 1998. (Also SUDMC Report 98--01.) 36. G.H. Golub and C.F. Van Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore MD, 3rd edition, 1996.
Pseudo-rigid and rigid-flexible bodies :,ii::~,iiiii''~i:,','i:,,iiii::i:,8i~ :::~:~' il::'~:::~:':iii!i!!'/' ii,i:::,i ,~,::i::ii'~:,'iii~, :i'~':i!''~::',i:~':i~',~"~i'~~,i~:',ii::~ii~ii ,i!'~,~,~ii~,:,~i,di'~i,:i'i~,:~:,:,i~~'i,:'~i,:,':ii,:i!i,Jii~i:~::il~~,'~,~i~,:i,'~!:,i',',i',i'i:,:~,i'::,':,:'i:,,i~:::iii i: i~'~!i:,:~:~':ii''i,i',!!',:,i:,, ,~i:'~::': ',':',i:,i'!,'~:~':,iii~:'~i'iii::i , :iii::::~i i~,i'':,:~i'i',:,::,~:~::i:~', ,i':!::i' ::::i'~':!::~::::~,i',i:,':~i:,~:':,~:::,:i:,~'~:, ,',~,ii','~::,'i,:,i~ii:,: i:i,::,::i:,:iii::~ii i ',:,i':,:,::'i~;,i:ii':,::,':i':',~::i:,:,::,~',i'i':i!'!i:,, !iiii', i:, :,~,il:,::::i:: : i ::i::ii: :,~ili::iii',i~'~~ii:,':,'~:, ::,iiii:,i',:,i!',:,i:,,~:,,~:~:,i:~, : ':'::,:'~','i:::i::ii i~:i::i:,ii::::~!i:ii:,ii',:,ii::,:,~,',',:i''~':,i::,!i::::i' '::,:~'::,iii:~:iil','i,:,!i'i':i',',i:il~:,',:!i,i':::,~i':i'',:if'i ::',i~:iil , :,i:,i:i::iiiiiiii:i:i:,,~:,'iliii:, '~i',:~,!:,i::':: :::ii::if:i~,::::i'ii:,i:'i::~::i:,!iilli'': i:::i'~ ,!:,::::,:~::i::ii' 'i: ::i'::i~iil
Many situations are encountered where treatment of the entire system as deformable bodies is neither necessary nor practical. For example, the frontal impact of a vehicle against a barrier requires a detailed modelling of the front part of the vehicle but the primary function of the engine and the rear part is to provide inertia, deformation being negligible for purposes of modelling the frontal impact. A second example, from geotechnical engineering, is the modelling of rock mass landslides or interaction between rocks on a conveyor belt where deformation of individual blocks is secondary. In this chapter we consider briefly the study of the first class of problems and in the next chapter the second type in much more detail. The above problem classes divide themselves into two further subclasses: one where it is necessary to include some simple mechanisms of deformation in each body (e.g. an individual rock piece) and the second in which the individual bodies have no deformation at all. The first class is calledpseudo-rigid body deformationl and the second rigid body behaviour. ~-Here we wish to illustrate how such behaviour can be described and combined in a finite element system. For the modelling of pseudo-rigid body analyses we follow closely the work of Cohen and Muncaster ~ and the numerical implementation proposed by Solberg and Papadopoulos. 3 The literature on rigid body analysis is extensive, and here we refer the reader to papers for additional details on methods and formulations beyond those covered here. 4-21
!iiiiiiiiiiiii~iiiii~iiiiiiiiiii~i~iiiiiiiiiiii~iiii~iii!!~i~iiiiiiiiiii~iiiiiii~iiii!ii!ii~iii
iiiiiiiiiiiiiiiiiiiii!~iii~iiiiiiiiiiiiiiiiiiiii!i!iiiiiiiiiiiiiiii!iiiiii!!iiiiiiiiiiiiiii i i iiiiii!i!iiiiiii i i iiiiiiiiiiiiiiiiii iiii iiiiiiiii!i iiii i!iiiiiiiii!ii!ii i !iiiiiii!iiiiii iiiiiiiiiiiiiiiiiiiiiiii i !ii
~T~,,~.~ v.,~L~',~~,,,,seuao-rlgm~motlOnS~'~,~
','~',~,~',~',~,~ ~~
~',~
;~ ~', ~ ~'~',~ ~ ~',~,~',~,~,~,~ ~ ~'~', ,~',
In this section we consider the analysis of systems which are composed of many small bodies, each of which is assumed to undergo large displacements and a uniform deformation.* The individual bodies which we consider are of the types shown in Fig. 8.1. In particular, a faceted shape can be constructed directly from a finite element discretization in which the elements are designated as all belonging to a single solid object or the individual bodies can be described by simple geometric forms such as discs or ellipsoids. * Higher-order body.
approximations
can
be
included
using
polynomial
approximation
for
the
deformation
of each
Pseudo-rigid motions 229
(a)
(b)
Fig. 8.1 Shapesfor pseudo-rigid and rigid body analysis: (a) ellipsoid; (b) faceted body. A homogeneous motion of a body may be written as
r
t) = ri(t) -+- Fii(t) [XI - RI]
(8.1)
in which XI is position, t is time, RI is some reference point in the undeformed body, ri is the position of the same point in the deformed body, and Fil is a constant deformation gradient. We note immediately that at time zero the deformation gradient is the identity tensor (matrix) and Eq. (8.1) becomes
r
O) - -
ri(O) -+- tSiI [ X I -- RI] = ri(O) + t ~ i l X l -- tSilRl ~ t S i l X l
(8.2)
where ri(O) "-- t~iiRi by definition. The behaviour of solids which obey the above description is sometimes referred to as analysis ofpseudo-rigid bodies. 1 A treatment by finite elements has been considered by Solberg and Papadopoulos, 3 and an alternative expression for motions restricted to incrementally linear behaviour has been developed by Shi, and the method is commonly called discontinuous deformation analysis (DDA). 22 The DDA form, while widely used in the geotechnical community, is usually combined with a simple linear elastic constitutive model and linear strain-displacement forms which can lead to large errors when finite rotations are encountered. Once the deformation gradient is computed, the procedures for analysis follow the methods described in Chapter 5. It is, of course, necessary to include the inertial term for each body in the analysis. No difficulties are encountered once a shape of each body is described and a constitutive model is introduced. For elastic behaviour it is not necessary to use a complicated model, and here use of the St Venant-Kirchhoff relation is adequate - indeed, if large deformations occur within an individual body the approximation of homogeneous deformation generally is not adequate to describe the solution. The primary difficulty for this class of problems is modelling the large number of interactions between bodies by contact phenomena and here the reader is referred to Chapters 7 and 9 and references on the subject for additional information on contact and other details. 22'23
230
Pseudo-rigidand rigid-flexible bodies
The pseudo-rigid body form can be directly extended to rigid bodies by using the polar decomposition on the deformation tensor. The polar decomposition of the deformation gradient may be given as 24-26 Fil -- A i j U j I
where
AilAij
"- (~lJ
and
AilAjl
= ~Sij
(8.3)
Here AiI is a rigid rotation* and UIj is a stretch tensor (which has eigenvalues Am as defined in Chapter 5). In the case of rigid motions the stretches are all unity and UIj simply becomes an identity. Thus, a rigid body motion may be specified as (gi(Xl, t) = ri(t) + A i l ( t ) [XI - RI]
(8.4)
q~(X, t) = r(t) + A ( t ) [ X - R]
(8.5)
or, in matrix form, as Alternatively, we can express the rigid motion using Eq. (8.1) and impose constraints to make the stretches unity. For example, in two dimensions we can represent the motion in terms of the displacements of the vertices of a triangle and apply constraints that the lengths of the triangle sides are unchanged during deformation. The constraints may be added as Lagrange multipliers or other constraint methods and the analysis may proceed directly from a standard finite element representation of the triangle. Such an approach has been used in reference 27 with a penalty method used to impose the constraints. Here we do not pursue this approach further and instead consider direct use of rigid body motions to construct the formulation. For subsequent use we note the form of the variation of a rigid motion and its incremental part. These may be expressed as A
3~b = 3r + 30A [X - R] dO - dr + dOA [X - R]
Using Eq. (8.5) these may be simplified to 3q~ - 3 r - ~30 dc~ - d r - ~dO
where
y = x- r (8.6)
where dO and ~0 are incremental and variational rotation vectors, respectively. In a similar manner we obtain the velocity for the rigid motion as -/"-
~w
(8.7)
in which f is translational velocity and w angular velocity, both at the centre of mass. The angular velocity is obtained by solving h = a,A * Often literature denotes this rotation as Rii; however, here we use to avoid confusion use Ail to denote rotation.
(8.8) RI
as a position of a point in the body and
Rigid motions 231 or
A-
A~
(8.9)
where f~ is the reference configuration angular velocity. 8 This is clearer by writing the equations in indicial form given by ~kil -- w i j A j l
(8.10)
-- A i j ~ ' 2 j l
where the velocity matrices are defined in terms of vector components and give the skew symmetric form Wij =
W3 --W 2
0
--W1
W1
0
(8.11)
and similarly for f2lj. The above form allows for the use of either the material angular velocity or the spatial one. Transformation between the two is easily performed since the rigid rotation must satisfy the orthogonality conditions ATA = A A T = I
(8.12)
at all times. Using Eqs (8.8) and (8.9) we obtain -- A ~ A T
(8.13)
= ATcDA
(8.14)
or by transforming in the opposite way
8.3.1 Equations of motion for a rigid body If we consider a single rigid body subjected to concentrated loads fa applied at points whose current position is Xa and locate the reference position for R at the centre of mass, the equations of equilibrium are given by conservation of linear momentum = ~
fa = f;
p = mr
(8.15)
a
where p defines a linear momentum, f is a resultant force and total mass of the body is computed from P
m - ]
J.
P0 dV
(8.16)
and conservation of angular momentum ~" - - Z
(Xa - -
r) • fa -- m;
Jr = I[w
(8.17)
a
where rr is the angular momentum of the rigid body, m is a resultant couple and ]I is the spatial inertia tensor.
232
Pseudo-rigidand rigid-flexible bodies The spatial inertia tensor (matrix) ITis computed from ]I = A J A T
(8.18)
where J is the inertia tensor (matrix) computed from an integral on the reference configuration and is given by ,~ - J~ P0 [(YTY) I -- yyT] dV
where
Y - X - R
(8.19)
Thus, description of an individual rigid body requires locating the centre of mass R and computing the total mass m and inertia matrix Jl. It is then necessary to integrate the equilibrium equations to define the position r and the orientation of the body A.
8.3.2 Construction from a finite element model If we model a body by finite elements, as described throughout this volume, we can define individual bodies or parts of bodies as being rigid. For each such body (or part of a body) it is then necessary to define the total mass, inertia matrix, and location of the centre of mass. This may be accomplished by computing the integrals given by Eqs (8.16) and (8.19) together with the relation to determine the centre of mass given by mR -- f n poX d V
(8.20)
In these expressions it is necessary only to define each point in the volume of an element by its reference position interpolation X. For solid (e.g. brick or tetrahedral) elements such interpolation is given by Eq. (5.50a) which in matrix form becomes (omitting the summation symbol) X = N.X~
(8.21)
This interpolation may be used to determine the volume element necessary to carry out all the integrals numerically by quadrature. 28 The total mass may now be computed as
m=~-'~(f~ podV) e
(8.22)
e
where ~"2 e is the reference volume of each element e. Use of Eq. (8.21) in Eq. (8.20) to determine the centre of mass now gives
R - --I _ m
e
PoN,~ dV) X~
(8.23)
e
and finally the reference inertia tensor (matrix) as J = E
M~;~ [(YSY~)I- Y~Y~] ; e
Y~ = X~ - R
(8.24)
Rigid motions 233 where
Mea3 = f~ PoN"N~dV
(8.25)
e
The above definition of Y, tacitly assumes that y ~ N, - 1. If other interpolations are used to define the shape functions (e.g. hierarchical shape functions) it is necessary to modify the above procedure to determine the mass and inertia matrix.
8.3.3 Transient solutions The integration of the translational rigid term r may be performed using any of the methods described in reference 28 or indeed by other methods described in the literature. The integration of the rotational part can also be performed by many schemes; however, it is important that updates of the rotation produce discrete time values for rigid rotations which retain an orthonormal character, that is, the An must satisfy the orthogonality condition given by Eq. (8.12). One procedure to obtain this is to assume that the angular velocity within a time increment is constant, being measured as
w(t) ~
(,~t3n+c~ =
•At
(8.26)
in which At is the time increment between tn and tn+l, 0 is the increment of rotation during the time step, and 0 < c~ < 1. The approximation (at,~nq_c~ --
(8.27)
(1 -- OL)(,~3 n -'1- OL(aL3n+ 1
is used to define intermediate values in terms of those at tn and tn+l. Equation (8.8) now becomes a constant coefficient ordinary differential equation which may be integrated exactly, yielding the solution A(t) -- exp[0(t - tn)/At]An
tn < t < tn-t-1
(8.28)
In particular at tn+~ we obtain An+c~ - exp[ce0]An This may also be performed using the material angular velocity f~.8 Many algorithms exist to construct the exponential of a matrix, and the closed-form expression given by the classical formula of Euler and Rodrigues (e.g. see Wittaker 29) is quite popular. This is given by
exp[0] - I +
1010 + _1 sin 2 101/202 101 2 [1012/2]
sin
where
101 =
[0T0 ] 1/2
(8.29)
This update may also be given in terms of quaternions and has been used for integration of both rigid body motions as well as for the integration of the rotations appearing in three-dimensional beam formulations (see Chapter 17). 8'30'31Another alternative to the
234
Pseudo-rigidand rigid-flexible bodies direct use of the exponential update is to use the Cayley transform to perform updates for A which remain orthonormal. Once the form for the update of the rigid rotation is defined any of the integration procedures defined in reference 28 may be used to advance the incremental rotation by noting that 0 or 19 (the material counterpart) are in fact the change from time t~ to tn+~. The reader also is referred to reference 8 for additional algorithms directly based on the GN 11 and GN22 methods. 28 Here forms for conservation of linear and angular momentum are of particular importance.
In some analyses the rigid body is directly attached to flexible body parts of the problem [Fig. 8.2(a)]. Consider a rigid body that occupies the part of the domain denoted a s ~2 r and is 'bonded' to a flexible body with domain f2f. In such a case the formulation to 'bond' the surface may be performed in a concise manner using Lagrange multiplier constraints. We shall find that these multiplier constraints can be easily eliminated from the analysis by a local solution process, as opposed to the need to carry them to the global solution arrays as was the case in their use in contact problems (see Sec. 5.3).
8.4.1 Lagrange multiplier constraints A simple two-dimensional rigid-flexible body problem is shown in Fig. 8.2(a) in which the interface will involve only three-nodal points. In Fig. 8.2(b) we show an exploded view between the rigid body and one of the elements which lies along the rigid-flexible interface. Here we need to enforce that the position of the two interface nodes for the element will have the same deformed position as the corresponding point on the rigid body. Such a constraint can easily be written using Eq. (8.4) as Ca = r(t) + A(t) [X~ - R] - x~(t) = 0
(8.30)
in which the subscript c~ denotes a node number. We can now modify a functional to include the constraint using a classical Lagrange multiplier approach in which we add the term 1-Irf = )~C~ = ,,~ [x~(t) - r(t) - A(t) [X~ - R]]
(8.31)
0
(a)
(b)
Fig. 8.2 Lagrange multiplier constraint between flexible and rigid bodies: (a)rigid-flexible body; (b) Lagrange multipliers.
Connecting a rigid body to a flexible body 235 Taking the variation we obtain
aFlrf = aA~ [x~ - r - A [X~ - R]] + A~ [ax~ - 6r - a0A [X~ - R]]
(8.32)
From this we immediately obtain the constraint equation and a modification to the equilibrium equations for each flexible node and the rigid body. Accordingly, the modified variational principle may now be written for a typical node a on the interface of the rigid body as ~I'I "4- ~I'Irf :
[~Xtz
6x~
~r
dO ~X,]
M.~v~ + M ~ v ~ + P. - f~, M ~ v ~ + M~av/~ + P~ - f~ + X~ (8.33)
#-f-X~
=0
e-m-yTA x~ - r -
A [ X ~ - R]
where y~ = x~ - r are the nodal values of y, fi are any other rigid body nodes connected to node a and #, v are flexible nodes connected to node a. Since the parameters x~ enter the equations in a linear manner we can use the constraint equation to eliminate their appearance in the equations. Accordingly, from the variation of the constraint equation we may write
6x~ = [I, y~]
(8.34)
60
which permits the remaining equations in Eq. (8.33) to be rewritten as
Mv~v~ + Mv~'b~ + P~ - f~ ~1"1 + ~I-Irf -- [~Xtt
&
~0]
p -- f + M~jr
+ M~av~ + P~ - f~
= 0
"/r - m - ~,T(M,~,% + M,~zi'a + P~ - f~)
(8.35) For use in a Newton solution scheme it is necessary to linearize Eq. (8.35). This is easily achieved
d(~I'I)+d(~I-Irf)= [ax~
ar I
30]
(K~)T
(Ktzfl)T
~X~
0
(Kav)T (Ka/~)T K.~ --yT(Kav)T --~T(Ka/~)T 0
qx~ KTo
dr
dO
(8.36)
236 Pseudo-rigidand rigid-flexible bodies Once again this form may be reduced using the equivalent of Eq. (8.34) for an incremental dx~ to obtain
d(~Fl) + d(~l-Irf )
[~x,
8r ~0]
r
(K~)T --(Kvr162 (KIzv)T - (K,~)T~'~ / (K~v)T K.~ + (K~3)T ,,T K [K~ + L--YT(K~.)T -Y~( ~;~)T
x x
l |
(8.37)
dr
dO Combining all the steps we obtain the set of equations for each rigid body as I
-- (Ktzfl)T~r/3
(Kvv)T
(Kv~)T
(K~)T
[K~ + (K~)T]
_~,T (Kav)T
_~,T (Ka~)T
]fdxv 1 ~d;
--(Ka/~)T~r/3 / [K~ + :9~(K~)T:9~IJ
(8.38)
=
m - ~ + ~'T'I'~ OL
in which ~ and ~ . are the residuals from the finite element calculation at node a and #, respectively. We recall from Chapter 5 that each is given by a form (8.39) which is now not zero since total balance of momentum includes the addition of the
~a. The above steps to compute the residual and the tangent can be performed in each element separately by noting that
,,ks = ~
Xe
(8.40)
e
where Xe denotes the contribution from element e. Thus, the steps to constrain a flexible body to a rigid body are once again a standard finite element assembly process and may easily be incorporated into a solution system. The above discussion has considered the connection between a rigid body and a body which is modelled using solid finite elements (e.g. quadrilateral and hexahedral elements in two and three dimensions, respectively). It is also possible directly to connect beam elements which have nodal parameters of translation and rotation. This is easily performed if the rotation parameters of the beam are also defined in terms of the rigid rotation A. In this case one merely transforms the rotation to be defined relative to the reference description of the rigid body rotation and assembles the result directly into the rotation terms of the rigid body. If one uses a rotation for both the beam and the rigid body which is defined in terms of the global Cartesian reference configuration no transformation is required. Shells can be similarly treated; however, it is best then to define the shell directly in terms of three rotation parameters instead of only two at points where connection is to be performed. 32'33
Multibody coupling by joints
Often it is desirable to have two (or more) rigid bodies connected in some specified manner. For example, in Fig. 8.3 we show a disc connected to an arm. Both are treated as rigid bodies but it is desired to have the disc connected to the arm in such a way that it can rotate freely about the axis normal to the page. This type of motion is characteristic of many rotating machine connections and it as well as many other types of connections are encountered in the study of rigid body motions. 4'34 This type of interconnection is commonly referred to as a joint. In quite general terms joints may be constructed by a combination of two types of simple constraints: translational constraints and rotational constraints.
8.5.1 Translation constraints The simplest type of joint is a spherical connection in which one body may freely rotate around the other but relative translation is prevented. Such a situation is shown in Fig. 8.3 where it is evident the spinning disc must stay attached to the rigid arm at its axle. Thus it may not translate relative to the arm in any direction (additional constraints are necessary to ensure it rotates only about the one axis - these are discussed in Sec. 8.5.2). If a full translation constraint is imposed a simple relation may be introduced as C j - - x (a) - x (b) - - 0 (8.41) where a and b denote two rigid bodies. Thus, addition of the Lagrange multiplier constraint (8.42) l - [ j - - ,,~jT[x(a) -- X(b)] imposes the spherical joint condition. It is necessary only to define the location for the spherical joint in the reference configuration. Denoting this as Xj (which is common
Fig. 8.3 Spinning disc constrained by a joint to a rigid arm.
237
238
Pseudo-rigidand rigid-flexible bodies to the two bodies) and introducing the rigid motion yields a constraint in terms of the rigid body positions as l-lj -- )kT[r (a) + A ( a ) ( x j
J
- R (a)) - r (b) - A(b)(xj -- R(b))]
(8.43)
The variation and subsequent linearization of this relation yields the contribution to the residual and tangent matrix for each body, respectively. This is easily performed using relations given above and is left as an exercise for the reader. If the translation constraint is restricted to be in one direction with respect to, say, body a it is necessary to track this direction and write the constraint accordingly. To accomplish this the specific direction of the body a in the reference configuration is required. This may be computed by defining two points in space X l and X2 from which a unit vector V is defined by V=
X2-Xl
IXz - Xll
(8.44)
The direction of this vector in the current configuration, v, may be obtained using the rigid rotation for body a v = A(a)v (8.45) A constraint can now be introduced into the variational problem as
I-lj =/kj{VT(A(a))T[r (a) + A(a)(xj - R (a)) - r O) - A~
- R(b))]}
(8.46)
where, owing to the fact there is only a single constraint direction, the Lagrange multiplier is a scalar/~j and, again, Xj denotes the reference position where the constraint is imposed. The above constraints may also be imposed by using a penalty function. The most direct form is to perturb each Lagrange multiplier form by a penalty tenn. Accordingly, for each constraint we write the variational problem as 1
2kj AJ
~
1-Ij -- )~jCj
2
(8.47)
where it is immediately obvious that the limit kj ---> oo yields exact satisfaction of the constraint. Use of a large kj and variation with respect to )~j give
E 1 ]
~%j C j - ~ ) ~ j
=0
(8.48)
and may easily be solved for the Lagrange multiplier as
)~j = kjCj
(8.49)
which when substituted back into Eq. (8.47) gives the classical form
nj -
kj
[cj
]2
(8.50)
The reader will recognize that Eq. (8.47) is a mixed problem, whereas Eq. (8.50) is irreducible. An augmented Lagrangian form is also possible following the procedures described in reference 28 and used in Chapter 7 for contact problems.
Multibody coupling by joints 239
8.5.2 Rotation constraints A second kind of constraint that needs to be considered relates to rotations. We have already observed in Fig. 8.3 that the disc is free to rotate around only one axis. Accordingly, constraints must be imposed which limit this type of motion. This may be accomplished by constructing an orthogonal set of unit vectors V I in the reference configuration and tracking the orientation of the deformed set of axes for each body as V~ c) --
~ilA(C)Vl
for c = a, b v~Vj - 5ij
(8.51)
A rotational constraint which imposes that axis i of body a remains perpendicular to axis j of body b may then be written as (v!a))Tv t
(b) - j
vT(A(a))TA(b)vj
-- 0
(8.52)
Example 8.1: Revolute joint As an example, consider the situation shown for the disc in Fig. 8.3 and define the axis of rotation in the reference configuration by the Cartesian unit vectors EI (i.e. V I = EI). If we let the disc be body a and the arm body b the set of constraints can be written as (where v3 is axis of rotation) X (a) __ X (c)
Cj =
~,*l[w(a)~Tw(b)J "3 = 0
(8.53)
l'w(a)~Tw(b) k*2 I *3
and included in a formulation using a Lagrange multiplier form
rlj = )~jT Cj
(8.54)
The modifications to the finite element equations are obtained by appending the variation and linearization of Eq. (8.54) to the usual equilibrium equations. Here five Lagrange multipliers are involved to impose the three translational constraints (spherical joint) and the angle constraints for the rotating disc. The set of constraints is known as a revolute joint, e
8.5.3 Library of joints Translational and rotational constraints may be combined in many forms to develop different types of constraints between rigid bodies. For the development it is necessary to have only the three types of constraints described above. Namely, the spherical joint, a single translational constraint, and a single rotational constraint. Once these are available it is possible to combine them to form classical constraint joints and here the reader is referred to the literature for the many kinds commonly encountered. 2,4,7,35 The only situation that requires special mention is the case when a series of rigid bodies is connected together to form a closed loop. In this case the method given above can lead to situations in which some of the joints are redundant. Using Lagrange multipliers this implies the resulting tangent matrix will be singular and, thus, one
240
Pseudo-rigid and rigid-flexible bodies cannot obtain solutions. Here a penalty method provides a viable method to circumvent this problem. The penalty method introduces elastic deformation in the joints and in this way removes the singular problem. If necessary an augmented Lagrangian method can be used to keep the deformation in the joint within required small tolerances. An alternative to this is to extract the closed loop rigid equations from the problem and use singular valued decomposition 36 to identify the redundant equations. These may then be removed by constructing a pseudo-inverse for the tangent matrix of the closed loop. This method has been used successfully by Chen to solve single loop problems. 35
'iii!iiiiiii iiiiiiiiiiii'ii,i'~i~ ,i~:~:',~::~::~iiiiiii'i,~:~i : :::~~i~,i,i:~,':ii':~iiii'i',i{iiii::~::::~'~'i~ii'~,iiiii::::~i~i'i~il~,ii,i~','i~,'i',i~',il!ili'~,i',iliiiiiii! i i:,i:,iiiiiiii~,~,',i~i{i'i,{!'!,!~,i ~!',iiiiiiiiii 'i,i ~""i~',~{'~,'!i~~:~ !i'i:~i',~'i~i'i i'~,i'~i,iiii'i!i'i,'i ':,,i'i,~i~,'il~,i!iiiiliiiiiiii' i i~i~,i',i!'i,:'iii'',',:iii!iii~'~'':iiili',iii',ii!'iii'~i',i',!i~,:,il,i',:{:i''i,~,i,i:ii!:i'',:,iii",','i, iii'!~,:',:ii',ii'"i~,,'~ili:,l'i,'ii'i:{~iiiiii'','i:,,',i!!,i~,i:i~:,'i:':,i:ii:i:,:i:,i::iif'~,i,'!',!i!'ii~'i,i'::,'i:~',,i'::',!,!iii!ii{i:i:,i!iii:,!:,i!':!:i{::i','i:i"~,ii'iiii!':i'~,'~,:!,i~:,i':i',i:i',,i','!'~,,i'i:!i~,,',',ii:: i"i'':,~,i i'!i,'','!~ii',i:,i' : ii',i'::,{:iiii'i':,,',i!:,:,i:,ii:!ii:i:' '~!~,i::ii':i':~ !i~i:'i :i,i','':iii',i,,'~~!: ':ii': i'
8.6.1 Rotating disc As a first example we consider a problem for the rotating disc on a rigid arm which is attached to a deformable base as shown in Fig. 8.4. The finite element model is constructed from four-node displacement elements in which a St Venant-Kirchhoff material model is used for the elastic part. The elastic properties in the model are E = 10 000 and z~ = 0.25, with a uniform mass density P0 = 5 throughout. The disc and arm are made rigid by using the procedures described in this chapter. The disc is attached to the arm by means of a revolute joint with the constraints imposed using the Lagrange multiplier method. The rigid arm is constrained to the elastic support by using the local Lagrange multiplier method described in Sec. 8.4. The problem is excited by a constant vertical load applied at the revolute joint and a torque applied to spin the disc. Each load is applied for the first 10 units of time. The mesh and configuration are shown in Fig. 8.4(a). Deformed positions of the model are shown at 2.5 unit intervals of time in Fig. 8.4(b)-(h). A marker element shows the position of the rotating disc. The displacements at the revolute joint and the radial exterior point at the marker element location are shown in Fig. 8.5.
8.6.2 Beam with attached mass ....................................................................................
--: .....................................................................................................................................................................................
As a second example we consider an elastic cantilever beam with an attached end mass of rectangular shape. The beam is excited by a horizontal load applied at the top as a triangular pulse for two units of time. The rigid mass is attached to the top of the beam by using the Lagrange multiplier method described in Sec. 8.4 and here it is necessary to constrain both the translation and the rotation parameters of the beam. The beam is three dimensional and has an elastic modulus of E = 100 000 and a moment of inertia in both directions of I1~ = 122 = 12. The beam mass density is low, with a value of /90 = 0.02. The tip mass is a cube with side lengths 4 units and mass density/90 = 1. The shape of the beam at several instants of time is shown in Fig. 8.6 and it is clear that large translation and rotation is occurring and also that the rigid block is correctly following a constrained rigid body motion.
Numerical examples 241
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Fig. 8.4 Rigid-flexible model for spinning disc: (a) problem definition, solutions at time; (b) t = 2.5 units; (c) t - 5.0 units; (d) t - 7.5 units; (e) t -- 10.0 units; (f) t - 12.5 units; (g) t - 15.0 units; (h) t - 17.5 units.
242
Pseudo-rigid and rigid-flexible bodies
15, -,r
0.5
E o
0
r 03
-0.5
A
A
l!
h
1 I I---Horizontall 10~"'Vertical I
A
".......",..-'I .- ::, ,,:I,- 1 0
4_~,,-
L
---Vertical J _ . -1 10 20 0 Time (a)
-15 -20 30
40
/""%
-.
'
0
10
'
20 Time (b)
30
40
Fig. 8.5 Displacements for rigid-flexible model for spinning disc. Displacement at: (a) revolute; (b) disc rim.
(a)
(c)
(b)
(e)
(f)
(g)
(d)
(h)
(i)
Fig. 8.6 Cantilever with tip mass: (a) t -- 2 units; (b) t -- 4 units; (c) t -- 6 units; (d) t - 10 units; (e) t -- 12 units; (f) t - 14 units; (g) t - 16 units; (h) t - 18 units; (i) t - 20 units.
1. H. Cohen and R.G. Muncaster. The Theory of Pseudo-rigid Bodies. Springer, New York, 1988. 2. A.A. Shabana. Dynamics of Multibody Systems. John Wiley & Sons, New York, 1989. 3. J.M. Solberg and P. Papadopoulos. A simple finite element-based framework for the analysis of elastic pseudo-rigid bodies. International Journal for Numerical Methods in Engineering, 45:1297-1314, 1999. 4. D.J. Benson and J.O. Hallquist. A simple rigid body algorithm for structural dynamics programs. International Journal for Numerical Methods in Engineering, 22:723-749, 1986.
References 243 5. R.A. Wehage and E.J. Haug. Generalized coordinate partitioning for dimension reduction in analysis of constrained dynamic systems. Journal of Mechanical Design, 104:247-255, 1982. 6. A. Cardona and M. Geradin. Beam finite element nonlinear theory with finite rotations. International Journal for Numerical Methods in Engineering, 26:2403-2438, 1988. 7. A. Cardona, M. Geradin and D.B. Doan. Rigid and flexible joint modelling in multibody dynamics using finite elements. Computer Methods in Applied Mechanics and Engineering, 89:395418,91. 8. J.C. Simo and K. Wong. Unconditionally stable algorithms for rigid body dynamics that exactly conserve energy and momentum. International Journal for Numerical Methods in Engineering, 31:19-52, 1991. [Addendum: 33:1321-1323, (1992).] 9. H.T. Clark and D.S. Kang. Application of penalty constraints for multibody dynamics of large space structures. Advances in the Astronautical Sciences, 79:511-530, 1992. 10. G.M. Hulbert. Explicit momentum conserving algorithms for rigid body dynamics. Computers and Structures, 44:1291-1303, 1992. 11. M. Geradin, D.B. Doan and I. Klapka. MECANO: a finite element software for flexible multibody analysis. Vehicle System Dynamics, 22:87-90, 1993. Supplement issue. 12. S.N. Atluri and A. Cazzani. Rotations in computational solid mechanics. Archives of Computational Methods in Engineering, 2:49-138, 1995. 13. O.A. Bauchau, G. Damilano and N.J. Theron. Numerical integration of nonlinear elastic multibody systems. International Journal for Numerical Methods in Engineering, 38:2727-2751, 1995. 14. J.A.C. Ambr6sio. Dynamics of structures undergoing gross motion and nonlinear deformations: a multibody approach. Computers and Structures, 59:1001-1012, 1996. 15. R.L. Huston. Multibody dynamics since 1990. Applied Mechanics Reviews, 49:$35-$40, 1996. 16. O.A. Bauchau and N.J. Theron. Energy decaying scheme for non-linear beam models. Computer Methods in Applied Mechanics and Engineering, 134:37-56, 1996. 17. O.A. Bauchau and N.J. Theron. Energy decaying scheme for non-linear elastic multi-body systems. Computers and Structures, 59:317-331, 1996. 18. C. Bottasso and M. Borri. Energy preserving/decaying schemes for nonlinear beam dynamics using the helicoidal approximation. Computer Methods in Applied Mechanics and Engineering, 143:393-415, 1997. 19. O.A. Bauchau. Computational schemes for flexible, nonlinear multi-body systems. Multibody System Dynamics, 2:169-222, 1998. 20. O.A. Bauchau and C.L. Bottasso. On the design of energy preserving and decaying schemes for flexible nonlinear multi-body systems. Computer Methods in Applied Mechanics and Engineering, 169:61-79, 1999. 21. O.A. Bauchau and T. Joo. Computational schemes for non-linear elasto-dynamics. International Journal for Numerical Methods in Engineering, 45:693-719, 1999. 22. G.-H. Shi. Block System Modelling by Discontinuous Deformation Analysis. Computational Mechanics Publications, Southampton, 1993. 23. E.G. Petocz. Formulation and analysis of stable time-stepping algorithms for contact problems. PhD thesis, Department of Mechanical Engineering, Stanford University, Stanford, California, 1998. 24. M.E. Gurtin. An Introduction to Continuum Mechanics. Academic Press, New York, 1981. 25. L.E. Malvern. Introduction to the Mechanics of a Continuous Medium. Prentice-Hall, Englewood Cliffs, NJ, 1969. 26. J. Bonet and R.D. Wood. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press, Cambridge, 1997. ISBN 0-521-57272-X. 27. J.C. Garcia Orden and J.M. Goicolea. Dynamic analysis of rigid and deformable multibody systems with penalty methods and energy-momentum schemes. Computer Methods in Applied Mechanics and Engineering, 188:789-804, 2000.
244 Pseudo-rigidand rigid-flexible bodies 28. O.C. Zienkiewicz, R.L. Taylor and J.Z. Zhu. The Finite Element Method: Its Basis and Fundamentals. Butterworth-Heinemann, Oxford, 6th edition, 2005. 29. E.T. Wittaker. A Treatise on Analytical Dynamics. Dover Publications, New York, 1944. 30. J.H. Argyris and D.W. Scharpf. Finite elements in time and space. Nuclear Engineering and Design, 10:456--469, 1969. 31. A. Ibrahimbegovic and M. A1 Mikdad. Finite rotations in dynamics of beams and implicit timestepping schemes. International Journal for Numerical Methods in Engineering, 41:781-814, 1998. 32. J.C. Simo. On a stress resultant geometrically exact shell model. Part VII: Shell intersections with 5/6 DOF finite element formulations. Computer Methods in Applied Mechanics and Engineering, 108:319-339, 1993. 33. P. Betsch, F. Gruttmann and E. Stein. A 4-node finite shell element for the implementation of general hyperelastic 3d-elasticity at finite strains. Computer Methods in Applied Mechanics and Engineering, 130:57-79, 1996. 34. H. Goldstein. Classical Mechanics. Addison-Wesley, Reading, 2nd edition, 1980. 35. A.J. Chen. Energy-momentum conserving methods for three dimensional dynamic nonlinear multibody systems. PhD thesis, Department of Mechanical Engineering, Stanford University, Stanford, California, 1998. (Also SUDMC Report 98--01.) 36. G.H. Golub and C.F. Van Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore MD, 3rd edition, 1996.