Pseudo-Rigid and Rigid-Flexible Bodies

CHAPTER

Pseudo-Rigid and Rigid-Flexible Bodies

9

9.1 Introduction 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 modeling 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 modeling the frontal impact. A second example, from geotechnical engineering, is the modeling 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 this class of problems. 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 called pseudo-rigid body deformation [1] and the second rigid body behavior [2]. Here we wish to illustrate how such behavior can be described and combined in a finite element system. For the modeling of pseudo-rigid body analyses we follow closely the work of Cohen and Muncaster [1] 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].

9.2 Pseudo-rigid motions 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.1 The individual bodies which we consider are of the types shown in Fig. 9.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 1 Higher-order approximations can be included using polynomial approximation for the deformation of each body. Finite Elem Meth Solid Str Mech 7E, Volume 2. http://dx.doi.org/10.1016/B978-1-85617-634-7.00009-0 © 2014 Elsevier Ltd. All rights reserved.

277

278

CHAPTER 9 Pseudo-Rigid and Rigid-Flexible Bodies

(a)

(b)

FIGURE 9.1 Shapes for pseudo-rigid and rigid body analysis: (a) ellipsoid; (b) faceted body.

object or the individual bodies can be described by simple geometric forms such as disks or ellipsoids. A homogeneous motion of a body may be written as   φi (X I , t) = ri (t) + Fi I (t) X I − R I

(9.1)

in which X I is position, t is time, R I is some reference point in the undeformed body, ri is the position of the same point in the deformed body, and Fi I is a constant deformation gradient. We note immediately that at time zero the deformation gradient is the identity tensor (matrix) and Eq. (9.1) becomes   φi (X I , 0) = ri (0) + δi I X I − R I = ri (0) + δi I X I − δi I R I ≡ δi I X I (9.2) where ri (0) = δi I R I by definition. The behavior of solids which obey the above description is sometimes referred to as analysis of pseudo-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 behavior 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 behavior 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

9.3 Rigid motions

the solution. The primary difficulty for this class of problems is modeling the large number of interactions between bodies by contact phenomena such as those covered in Chapter 8. Here, however, many contacts are encountered and special methods can speed up the calculations, see for example, Refs. [22,23].

9.3 Rigid motions 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] Fi I = i J U J I

where

i I i J = δ I J

and

i I  j I = δi j

(9.3)

Here i I is a rigid and U I J is a stretch tensor (which has eigenvalues λm as defined in Chapter 5). In the case of rigid motions the stretches are all unity and U I J simply becomes an identity. Thus, a rigid body motion may be specified as   φi (X I , t) = ri (t) + i I (t) X I − R I (9.4) rotation2

or, in matrix form, as

  φ(X, t) = r(t) + (t) X − R

(9.5)

Alternatively, we can express the rigid motion using Eq. (9.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 Ref. [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 (see Chapter 11 for additional details). These may be expressed as    X − R δφ = δr + δθ    X − R dφ = dr + dθ where the matrix notation for the vector cross product introduced in Eqs. (5.124) and (5.125) is used. Using Eq. (9.5) these may be simplified to δφ = δr − yˆ δθ dφ = dr − yˆ dθ

(9.6)

where y = x − r and dφ and δθ are incremental and variational rotation vectors, respectively. 2 Often literature denotes this rotation as R ; however, here we use R as a position of a point in the iI I body and to avoid confusion use i I to denote rotation.

279

280

CHAPTER 9 Pseudo-Rigid and Rigid-Flexible Bodies

In a similar manner we obtain the velocity for the rigid motion as φ˙ = r˙ − yˆ ω

(9.7)

in which r˙ is translational velocity and ω angular velocity, both at the center of mass. The angular velocity is obtained by solving

or

˙ = ω ˆ 

(9.8)

˙ =  ˆ 

(9.9)

where  is the reference configuration angular velocity [8]. This is clearer by writing the equations in indicial form given by ˙ i I = ωi j  j I = i J  J I 

(9.10)

where the velocity matrices are defined in terms of vector components and give the skew-symmetric form [viz. Eq. (5.125)] ⎡ ⎤ 0 −ω3 ω2 ωi j = ⎣ ω3 0 −ω1 ⎦ (9.11) −ω2 ω1 0 and similarly for  I J . 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 T  = T = I

(9.12)

at all times. Using Eqs. (9.8) and (9.9) we obtain ˆ T ωˆ = 

(9.13)

or by transforming in the opposite way ˆ = T ω ˆ 

(9.14)

9.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 center of mass, the equations of equilibrium are given by conservation of linear momentum

fa = f, p = m r˙ (9.15) p˙ = a

where p defines a linear momentum, f is a resultant force, and the total mass of the body is computed from  ρ0 dV (9.16) m= 

9.3 Rigid motions

and conservation of angular momentum

π˙ = (xa − r) × fa = m,

π = Iω

(9.17)

a

where π is the angular momentum of the rigid body, m is a resultant couple, and I is the spatial inertia tensor. The spatial inertia tensor (matrix) I is computed from I = J T

(9.18)

where J is the inertia tensor (matrix) computed from an integral on the reference configuration and is given by 

J = where Y=X−R (9.19) ρ0 (YT Y)I − YYT dV 

Thus, description of an individual rigid body requires locating the center of mass R and computing the total mass m and inertia matrix J . It is then necessary to integrate the equilibrium equations to define the position r and the orientation of the body .

9.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 center of mass. This may be accomplished by computing the integrals given by Eqs. (9.16) and (9.19) together with the relation to determine the center of mass given by  ρ0 X dV (9.20) mR = 

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) ˜a X = Na X

(9.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 

 ρ0 dV (9.22) m= e

e

281

282

CHAPTER 9 Pseudo-Rigid and Rigid-Flexible Bodies

where e is the reference volume of each element e. Use of Eq. (9.21) in Eq. (9.20) to determine the center of mass now gives   1

R= ρ0 Na dV Xa (9.23) m e e and finally the reference inertia tensor (matrix) as  

J = Mab YaT Yb I − Ya YbT ,

Ya = Xa − R

(9.24)

e



where e Mab =

e

ρ0 Na Nb dV

(9.25)

 The above definition of Ya tacitly assumes that a Na = 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.

9.3.3 Transient solutions The integration of the translational rigid term r may be performed using any of the methods described in Ref. [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 n must satisfy the orthogonality condition given by Eq. (9.12). One procedure to obtain this is to assume that the angular velocity within a time increment is constant, being measured as ω(t) ≈ ωn+α =

1 θ

t

(9.26)

in which t is the time increment between tn and tn+1 , θ is the increment of rotation during the time step, and 0 ≤ α ≤ 1. The approximation ωn+α = (1 − α)ωn + αωn+1

(9.27)

is used to define intermediate values in terms of those at tn and tn+1 . Equation (9.8) now becomes a constant coefficient ordinary differential equation which may be integrated exactly, yielding the solution ˆ − tn )/ t]n (t) = exp[θ(t

tn ≤ t ≤ tn+1

(9.28)

In particular at tn+α we obtain n+α = exp[α θˆ]n This may also be performed using the material angular velocity  [8]. Many algorithms exist to construct the exponential of a matrix, and the closed-form expression

9.4 Connecting a rigid body to a flexible body

given by the classical formula of Euler and Rodrigues (e.g., see Whittaker [29]) is quite popular. This is given by ˆ =I+ exp[θ]

sin | θ | ˆ 1 sin2 | θ | /2 ˆ2 θ+ θ |θ | 2 [| θ |2 /2]

where



1/2 | θ |= θ T θ

(9.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 11) [8,30,31]. Another alternative to the direct use of the exponential update is to use the Cayley transform to perform updates for  which remain orthonormal. Once the form for the update of the rigid rotation is defined any of the integration procedures defined in Ref. [28] may be used to advance the incremental rotation by noting that θ or  (the material counterpart) is in fact the change from time tn to tn+1 . The reader also is referred to Ref. [8] for additional algorithms directly based on the generalized midpoint and Newmark methods. Here forms for conservation of linear and angular momentum are of particular importance.

9.4 Connecting a rigid body to a flexible body In some analyses the rigid body is directly attached to flexible body parts of the problem [Fig. 9.2a]. Consider a rigid body that occupies the part of the domain denoted as r and is “bonded” to a flexible body with domain  f . 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 Section 8.3).

9.4.1 Lagrange multiplier constraints A simple two-dimensional rigid-flexible body problem is shown in Fig. 9.2a in which the interface will involve only three-nodal points. In Fig. 9.2b we show an exploded

(a)

(b)

FIGURE 9.2 Lagrange multiplier constraint between flexible and rigid bodies: (a) rigid-flexible body; (b) Lagrange multipliers.

283

284

CHAPTER 9 Pseudo-Rigid and Rigid-Flexible Bodies

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 for a node xc using Eq. (9.5) as   (9.30) Cc = r(t) + (t) Xc − R − xc (t) = 0 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    (9.31) r f = λcT Cc = λcT xc (t) − r(t) − (t) Xc − R Taking the variation we obtain       δr f = δλcT xc − r −  Xc − R + λcT δxc − δr − δθ  Xc − R

(9.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 α on the interface of the rigid body as   δ + δr f = δxaT δxcT δr T δθ T δλcT ⎧ ⎫ Mab v˙ b + Mad v˙ d + Pa − fa ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ˙ ˙ v v M + M + P − f + λ cd d c c c⎪ ⎪ cb b ⎪ ⎨ ⎬ =0 (9.33) × p˙ − f − λc ⎪ ⎪ ⎪ ⎪ T ⎪ ⎪ ⎪ ⎪ π˙ − m − yˆ c λc ⎪ ⎪ ⎪ ⎪ ⎪ ⎪   ⎩ ⎭ xc − r −  Xc − R where yc = xc −r are the nodal values of y, d are any other rigid body nodes connected to node c, and a, b are all flexible nodes not connected to nodes c or d. Since the parameters xc 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   δr δxc = [I, yˆ c ] (9.34) δθ which permits the remaining equations in Eq. (9.33) to be rewritten as ⎧ ⎫ ⎪ Mab v˙ b + Mad v˙ d + Pa − fa ⎪ ⎬  T T T ⎨ p˙ − f + Mcb v˙ d + Mcd v˙ d + Pc − fc δ + δr f = δxa δr δθ =0 ⎪ ⎪ ⎩ ⎭ π˙ − m − yˆ cT (Mcb v˙ b + Mcd v˙ d + Pc − fc ) (9.35)

9.4 Connecting a rigid body to a flexible body

For use in a Newton solution scheme it is necessary to linearize Eq. (9.35). This is easily achieved:   d(δ) + d(δr f ) = δxaT δr T δθ T ⎫ ⎧ ⎡ ⎤ dxb ⎪ ⎪ ⎪ Kad 0 0 ⎪ Kab ⎪ ⎪ ⎬ ⎨ ⎢ ⎥ dxd p Kcd K 0 ⎦ (9.36) × ⎣ Kcb ⎪ ⎪ ⎪ ⎪ dr ⎪ ⎪ ⎭ −ˆycT Kcd −ˆycT Kcd 0 Kθ ⎩ dθ where Kab denotes a tangent stiffness term for the nodal pairs a, b of a flexible element. Once again this form may be reduced using the equivalent of Eq. (9.34) for an incremental dxd to obtain   d(δ) + d(δr f ) = δxaT δr T δθ T ⎫ ⎤⎧ ⎡ Kad −Kad yˆ d Kab dx ⎪ ⎬ ⎨ b⎪ ⎥ ⎢ K p + Kcd −Kcd yˆ d × ⎣ Kcb (9.37) ⎦ dr ⎪  θ  ⎪ ⎭ ⎩ T T T dθ K + yˆ c Kcd yˆ d −ˆyc Kcb −ˆyc Kcd Combining all the steps we obtain the set of equations for each rigid body as ⎫ ⎤⎧ ⎡ Kad −Kad yˆ d Kab dxb ⎪ ⎪ ⎬ ⎨  p  ⎥ ⎢ K + Kcd −Kcd yˆ d ⎦ dr ⎣ Kcb ⎪  θ  ⎪ ⎭ ⎩ dθ K + yˆ cT Kcd yˆ d −ˆycT Kcb −ˆycT Kcd ⎧ ⎫

⎪ ⎪ ⎨ a ⎬ = f − p˙ + c (9.38) ⎪ ⎪ ⎩ ⎭ m − π˙ + yˆ cT c in which c and a are the residuals from the finite element calculation at node c and a, respectively. We recall from Chapter 5 that each is given by a form ψ a = fa − Pa (σ ) − Mab v˙ b − Mad v˙ d

(9.39)

which is now not zero since the total balance of momentum includes the addition of the λc . The above steps to compute the residual and the tangent can be performed in each element separately by noting that

λec (9.40) λc = e

where λec 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.

285

286

CHAPTER 9 Pseudo-Rigid and Rigid-Flexible Bodies

FIGURE 9.3 Spinning disk constrained by a joint to a rigid arm.

The above discussion has considered the connection between a rigid body and a body which is modeled using solid finite elements (e.g., quadrilateral and hexahedral elements in two and three dimensions, respectively). It is also possible directly to connect beam (or rod) 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  (viz. Chapter 13). 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 (such as those considered in Chapter 14); 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].

9.5 Multibody coupling by joints Often it is desirable to have two (or more) rigid bodies connected in some specified manner. For example, in Fig. 9.3 we show a disk connected to an arm. Both are treated as rigid bodies but it is desired to have the disk 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. Many other types of connections also 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: translation constraints and rotation constraints.

9.5 Multibody coupling by joints

9.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. 9.3 where it is evident the spinning disk 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 Section 9.5.2). If a full translation constraint is imposed a simple relation may be introduced as (9.41) C j = x(a) − x(b) = 0 where a and b denote two rigid bodies. Thus, addition of the Lagrange multiplier constraint (9.42)  j = λTj [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 X j (which is common to the two bodies) and introducing the rigid motion yields a constraint in terms of the rigid body positions as  j = λTj [r(a) + (a) (X j − R(a) ) − r(b) − (b) (X j − R(b) )]

(9.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 body a in the reference configuration is required. This may be computed by defining two points in space X1 and X2 from which a unit vector V is defined by V=

X2 − X1 | X2 − X1 |

(9.44)

The direction of this vector in the current configuration, v, may be obtained using the rigid rotation for body a: (9.45) v = (a) V A constraint can now be introduced into the variational problem as  j = λ j {VT ((a) )T [r(a) + (a) (X j − R(a) ) − r(b) − (b) (X j − R(b) )]} (9.46) where, owing to the fact there is only a single constraint direction, the Lagrange multiplier is a scalar λ j and, again, X j denotes the reference position where the constraint is imposed.

287

288

CHAPTER 9 Pseudo-Rigid and Rigid-Flexible Bodies

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 term. Accordingly, for each constraint we write the variational problem as 1 2 j = λjCj − λ (9.47) 2k j j where it is immediately obvious that the limit k j → ∞ yields exact satisfaction of the constraint. Use of a large k j and variation with respect to λ j give   1 (9.48) δλ j C j − λ j = 0 kj and may easily be solved for the Lagrange multiplier as λj = kjCj

(9.49)

which when substituted back into Eq. (9.47) gives the classical form kj (9.50)  j = [C j ]2 2 The reader will recognize that Eq. (9.47) is a mixed problem, whereas Eq. (9.50) is irreducible. An augmented Lagrangian form is also possible following the procedures described in Ref. [28] and used in Chapter 8 for contact problems.

9.5.2 Rotation constraints A second kind of constraint that needs to be considered relates to rotations. We have already observed in Fig. 9.3 that the disk 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 vi(c) = δi I (c) V I for c = a, b vTI V J = δ I J

(9.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 (a)

(b)

(vi )T v j = VTI ((a) )T (b) V J = 0 Example 9.1.

(9.52)

Revolute joint

As an example, consider the situation shown for the disk in Fig. 9.3 and define the axis of rotation in the reference configuration by the Cartesian unit vectors E I (i.e., V I = E I ). If we let the disk be body a and the arm body b the set of constraints can be written as (where v3 is the axis of rotation) ⎫ ⎧ (a) x − x(c) ⎪ ⎪ ⎪ ⎪ ⎬ ⎨ (a) T (b) =0 (9.53) C j = (v1 ) v3 ⎪ ⎪ ⎪ ⎭ ⎩ (a) T (b) ⎪ (v2 ) v3

9.6 Numerical examples

and included in a formulation using a Lagrange multiplier form:  j = λTj C j

(9.54)

The modifications to the finite element equations are obtained by appending the variation and linearization of Eq. (9.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 disk. The set of constraints is known as a revolute joint [2].

9.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 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].

9.6 Numerical examples 9.6.1 Rotating disk As a first example we consider a problem for the rotating disk on a rigid arm which is attached to a deformable base as shown in Fig. 9.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 ν = 0.25, with a uniform mass density ρ0 = 5 throughout. The disk and arm are made rigid by using the procedures described in this chapter. The disk 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

289

290

CHAPTER 9 Pseudo-Rigid and Rigid-Flexible Bodies

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

FIGURE 9.4 Rigid-flexible model for spinning disk: (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.

9.6 Numerical examples

(a)

(b)

FIGURE 9.5 Displacements for rigid-flexible model for spinning disk. Displacement at: (a) revolute; (b) disk rim.

by using the local Lagrange multiplier method described in Section 9.4. The problem is excited by a constant vertical load applied at the revolute joint and a torque applied to spin the disk. Each load is applied for the first 10 units of time. The mesh and configuration are shown in Fig. 9.4a. Deformed positions of the model are shown at 2.5 unit intervals of time in Fig. 9.4b–h. A marker element shows the position of the rotating disk. The displacements at the revolute joint and the radial exterior point at the marker element location are shown in Fig. 9.5.

9.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 Section 9.4 and here it is necessary to constrain both the translation and 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 I11 = I22 = 12. The beam mass density is low, with a value of ρ0 = 0.02. The tip mass is a cube with side lengths 4 units and mass density ρ0 = 1. The shape of the beam at several instants of time is shown in Fig. 9.6 and it is clear that large translation and rotation is occurring and that the rigid block is correctly following a constrained rigid body motion.

9.6.3 Biofidelic rear impact dummy The Abaqus BioRID II model shown in Fig. 9.7 is used to evaluate the risk of neck injury during a rear-end collision. The Abaqus model, along with the physical test dummy it represents, is used by automotive original equipment manufacturers and seat suppliers to design head restraint systems that minimize whiplash injuries that can often occur in low-speed rear collisions.

291

292

CHAPTER 9 Pseudo-Rigid and Rigid-Flexible Bodies

(a)

(e)

(b)

(c)

(f)

(g)

(d)

(h)

(i)

FIGURE 9.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.

FIGURE 9.7 Biofidelic rear impact dummy. Image courtesy of Dassault Systèmes SIMULIA.

9.7 Concluding remarks

FIGURE 9.8 Sorting of randomly sized particles. Image courtesy of Dassault Systèmes SIMULIA.

9.6.4 Sorting of randomly sized particles When many small rigid or pseudo-rigid bodies are involved one can use analyses based on large numbers of discrete elements. The discrete element method (DEM) [37–40] is commonly used for such analyses. A main feature of such analyses is rapid treatment of contact interactions. The image in Fig. 9.8 shows the application of DEM in Abaqus, where spherical particles with randomly varying radii are sorted by screens with progressively smaller openings.

9.7 Concluding remarks In this chapter we have considered the solution of problems in which a portion or all parts may be treated as a rigid body. In addition the possibility of including simple deformation modes is addressed. In principle this class of problems belongs to a class of reduced order models for which much current research is taking place. Combining rigid or pseudo-rigid parts with deformable treatment as covered in earlier chapters permits the solution of a vary wide class of problems as illustrated above by the example of the rear impact dummy. When many rigid or pseudo-rigid bodies are involved special analysis methods are needed. One is the discrete element method (DEM) cited in the last example. For cases where spherical objects or clusters of spherical objects are needed one can consider use of discrete deformation analysis (DDA) [22,41,42].

293

294

CHAPTER 9 Pseudo-Rigid and Rigid-Flexible Bodies

In the next chapters we consider an important class of problems in which deformable models have a reduction in dimension, but only in one or two directions. This is the structural class of beam, plate, and shell problems. In addition to considering a new class of problems we also introduce a more modern approach to computational mechanics—one that has foundations in both mechanics and mathematics. This is extremely important to gain an ability to follow much of the recent research related to computational methods.

References [1] H. Cohen, 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, P. Papadopoulos, A simple finite element-based framework for the analysis of elastic pseudo-rigid bodies, Int. J. Numer. Methods Eng. 45 (1999) 1297–1314. [4] D.J. Benson, J.O. Hallquist, A simple rigid body algorithm for structural dynamics programs, Int. J. Numer. Methods Eng. 22 (1986) 723–749. [5] R.A. Wehage, E.J. Haug, Generalized coordinate partitioning for dimension reduction in analysis of constrained dynamic systems, J. Mech. Des. 104 (1982) 247–255. [6] A. Cardona, M. Geradin, Beam finite element nonlinear theory with finite rotations, Int. J. Numer. Methods Eng. 26 (1988) 2403–2438. [7] A. Cardona, M. Geradin, D.B. Doan, Rigid and flexible joint modelling in multibody dynamics using finite elements, Comput. Methods Appl. Mech. Eng. 89 (1991) 395–418. [8] J.C. Simo, K. Wong, Unconditionally stable algorithms for rigid body dynamics that exactly conserve energy and momentum, Int. J. Numer. Methods Eng. 31 (1991) 19–52 (Addendum: 33 (1992) 1321–1323). [9] H.T. Clark, D.S. Kang, Application of penalty constraints for multibody dynamics of large space structures, Adv. Astronaut. Sci. 79 (1992) 511–530. [10] G.M. Hulbert, Explicit momentum conserving algorithms for rigid body dynamics, Comput. Struct. 44 (1992) 1291–1303. [11] M. Geradin, D.B. Doan, I. Klapka, MECANO: a finite element software for flexible multibody analysis, Veh. Syst. Dyn. 22 (Suppl.) (1993) 87–90. [12] S.N. Atluri, A. Cazzani, Rotations in computational solid mechanics, Arch. Comput. Methods Eng. 2 (1995) 49–138. [13] O.A. Bauchau, G. Damilano, N.J. Theron, Numerical integration of nonlinear elastic multibody systems, Int. J. Numer. Methods Eng. 38 (1995) 2727–2751. [14] J.A.C. Ambrósio, Dynamics of structures undergoing gross motion and nonlinear deformations: a multibody approach, Comput. Struct. 59 (1996) 1001–1012.

References

[15] R.L. Huston, Multibody dynamics since 1990, Appl. Mech. Rev. 49 (1996) S35–S40. [16] O.A. Bauchau, N.J. Theron, Energy decaying scheme for non-linear beam models, Comput. Methods Appl. Mech. Eng. 134 (1996) 37–56. [17] O.A. Bauchau, N.J. Theron, Energy decaying scheme for non-linear elastic multi-body systems, Comput. Struct. 59 (1996) 317–331. [18] C. Bottasso, M. Borri, Energy preserving/decaying schemes for nonlinear beam dynamics using the helicoidal approximation, Comput. Methods Appl. Mech. Eng. 143 (1997) 393–415. [19] O.A. Bauchau, Computational schemes for flexible, nonlinear multi-body systems, Multibody Syst. Dyn. 2 (1998) 169–222. [20] O.A. Bauchau, C.L. Bottasso, On the design of energy preserving and decaying schemes for flexible nonlinear multi-body systems, Comput. Methods Appl. Mech. Eng. 169 (1999) 61–79. [21] O.A. Bauchau, T. Joo, Computational schemes for non-linear elasto-dynamics, Int. J. Numer. Methods Eng. 45 (1999) 693–719. [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, Ph.D. Thesis, Department of Mechanical Engineering, Stanford University, Stanford, California, 1998. [24] M.E. Gurtin, An introduction to continuum mechanics, Mathematics in Science and Engineering, vol. 158, 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, R.D. Wood, Nonlinear Continuum Mechanics for Finite Element Analysis, Cambridge University Press, Cambridge, 1997, ISBN 0-521-57272-X. [27] J.C. García Orden, J.M. Goicolea, Dynamic analysis of rigid and deformable multibody systems with penalty methods and energy-momentum schemes, Comput. Methods Appl. Mech. Eng. 188 (2000) 789–804. [28] O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, seventh ed., Elsevier, Oxford, 2013. [29] E.T. Wittaker, A Treatise on Analytical Dynamics, Dover Publications, New York, 1944. [30] J.H. Argyris, D.W. Scharpf, Finite elements in time and space, Nucl. Eng. Des. 10 (1969) 456–469. [31] A. Ibrahimbegovic, M. Al Mikdad, Finite rotations in dynamics of beams and implicit time-stepping schemes, Int. J. Numer. Methods Eng. 41 (1998) 781–814. [32] J.C. Simo, On a stress resultant geometrically exact shell model. Part VII: Shell intersections with 5/6 DOF finite element formualtions, Comput. Methods Appl. Mech. Eng. 108 (1993) 319–339.

295

296

CHAPTER 9 Pseudo-Rigid and Rigid-Flexible Bodies

[33] P. Betsch, F. Gruttmann, E. Stein, A 4-node finite shell element for the implementation of general hyperelastic 3D-elasticity at finite strains, Comput. Methods Appl. Mech. Eng. 130 (1996) 57–79. [34] H. Goldstein, Classical Mechanics, second ed., Addison-Wesley, Reading, MA, 1980. [35] A.J. Chen, Energy-momentum conserving methods for three dimensional dynamic nonlinear multibody systems, Ph.D. Thesis, Department of Mechanical Engineering, Stanford University, Stanford, California, 1998 (also SUDMC Report 98–01). [36] G.H. Golub, C.F. Van Loan, Matrix Computations, third ed., The Johns Hopkins University Press, Baltimore, MD, 1996. [37] P.A. Cundall, D.L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1979) 47–65. [38] C. Thornton, Applications of DEM to process engineering problems, Eng. Comput. 9 (1992) 289–297. [39] P.W. Cleary, Large scale industrial DEM modelling, in: Discrete Element Methods: Numerical Modeling of Discontinua, 3rd International Conference on Discrete Element Methods, ASCE Geotechnical Special Publication No. 117, The Geo-Institute of the American Society of Civil Engineers, 2002. [40] O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, sixth ed., Elsevier, Oxford, 2005. [41] L. Jing, Formulation of discontinuous deformation (DDA) - an implicit discrete element model for block systems, Eng. Geol. 49 (1998) 371–381. [42] C.J. Pearce, A. Thavalingam, Z. Liao, N. Biˇcaniˇc, Computational aspects of the discontinuous deformation analysis framework for modelling concrete, Eng. Fract. Mech. 65 (2000) 283–298.