The implementation of a massless deformable body in MADYMO

The implementation of a massless deformable body in

MADYMO

Citation for published version (APA):

Vosbeek, P. H. J. (1990). The implementation of a massless deformable body in MADYMO. Technische Universiteit Eindhoven.

Document status and date: Published: 01/01/1990 Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers) Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

EINDHOVEN UNIVERSITY OF TECHNOLOGY

Department of Mathematics and Computing Science

The implementation of

a massless deformable body

inMADYMO

by

P.H.J. Vosbeek

Preface

After my graduation at the Eindhoven University of Technology, I was able to continue my investigations at the TNO Road-Vehicles Research Institute, Delft, for two months. This report is the result of these investigations. It describes the implementation of a massless deformable body in MADYMO, and gives two examples. The third and fourth chapter of my Master's thesis were rewritten and another chapter was added. The latter chapter describes the implementation of a linear-viscoelastic beam in MADYMO.

Februari 1990, Pieter Vosbeek.

Contents

Preface

Contents

1 Introduction

2 Two rigid bodies connected by a massless deformable body 2.1 Introduction

2.2 The relative displacement and velocity of the two rigid bodies 2.3 The relative orientation and angular velocity of the two

rigid bodies

2.4 The relative orientation in a linear theory

2.5 The forces and torques applied to the two rigd bodies 2.6 The implementation in MADYMO

2.7 Input description

3 A linear elastic beam with a circular cross-section 3.1 Introduction

3.2 The loads on the linear elastic beam 3.3 Some simulations

4 A linear visco-elastic beam with a circular cross-section 4.1 Introduction. Basic equations

1 2 4 6 6 6 10 12 14 17 21 24 24 24 27 32 32

4.2

Constitutive equations

33

4.3

Boundary conditions

34

4.4

The transformation of the problem

35

4.5

The solution of the problem

36

4.6

Some simple visco--€lastic models

38

4.7

The solution for the standard linear solid

42

4.8

The formulation in ODE's

44

4.9

The implementation in MADYMO

47

4.10

Some simulations

47

5 Conclusions and suggestions for further improvements 53

1 Introduction

MADYMO is a compact computer program for simulating the dynamic motion of systems of connected rigid bodies. The program has been developed especially for the simulation of the motion of the human body during an impact.

In MADYMO, the human body is represented by a number of connected rigid bodies. In order to describe the motion of some parts of the human body (like the neck and the back) accurately, these parts ought to be modelled by a large number of rigid bodies. However, this leads to a large system of equations of motion and it takes a considerable amount of time to solve such a system. An alternative is to represent such parts by deformable bodies. It is very complex to add the possibility of modelling deformable bodies to MADYMO. For reasons of simplicity the deformable body could be assumed to be massless. (Its mass could be taken into account by adding it to the masses of the rigid bodies the deformable body connects.) In the massless case the only contribution of a deformable body to. the equations of motion is a contribution to the right-hand side of the equations of motion (the forces and torques applied to the rigid bodies, due to the deformable body). These loads are functions of the relative position and orientation of the connected rigid bodies and their first time derivatives. They can be taken into account by writing a user-defined subroutine USERFO which adds the forces and torques to the right-hand side of the equations of motion.

In the next chapter the relative displacements and changes oforientation of two rigid bodies, connected by a massless deformable body, shall be described. Then, the loads the massless body exerts on the rigid bodies shall be determined. In the Chapters 3 and 4 two deformable bodies are presented: a linear elastic beam and a

linear visco-elastic beam. Both beams have a circular cross-section. The loads that these beams apply to the two rigid rigid bodies they connect are determined. In the last chapter of this report some improvements that have to be made in order to implement the developed subroutines in MADYMO shall be discussed.

2

Two rigid bodies connected by a massless deformable body

2.1 Introduction

Consider two rigid bodies connected by a massless deformable body. The deformable body is rigidly connected to each rigid body in one single point. In this chapter the. forces and torques the deformable body exerts on the two rigid bodies are determined. These loads are functions of the relative displacement and velocity of the points of connection and of the change of orientation and the relative angular velocity of the two rigid bodies. The latter quantities can be seen as the displacement, velocity, orientation and angular velocity of one rigid body, while the other is kept in a fixed position. In the next two sections expressions for these quantities are derived. In Section 2.4

a

vector is introduced that describes the relative orientation in a linear theory. By means of a relation between the relative displacement, velocity, orientation and angular velocity on one hand and the forces and torques exerted on the deformable body on the other, the loads applied to the two rigid bodies are calculated in Section 2.5. This relation depends on the material the massless body is made of and on its shape. The last two sections of this chapter describe the implementation of the theory in

MADYMO

and give an input description of the subroutine DEBOFO that calculates the forces and torques.

2.2 The relative displacement and velocity of the two rigid bodies

Denote the two rigid bodies by $ ₁ and $ ₂ and the massless deformable body by $. The body$ is connected to $ ₁ in the point .9'₁ and to .!'4₂ in .9'_2. The region G(t), occupied by the three bodies in the three-<limensional euclidian space IR₃ at the time

t,

is referred to as the configuration of .!fi _1, .!fi ₂ and .!fi. In order to describe the motion of the three bodies, a reference configuration G₀ is defined. This is a time-independent configuration in which all relevant quantities are known.

Choose in each of the rigid bodies .!fi i a fixed point .Ai with position vector Ai

with respect to the origin() of an inertial space, in the reference configuration G_0. Let.% E .!fi i be an arbitrary material point of the rigid body .!fi i and let X be its position in G_0• The position x( t) E G( t) of point.% at the time t can then be written

as,

(2.2.1)

Here, Ri ( t) is a proper orthogonal mapping (Rt ( t) Ri ( t)

=

Ri ( t)

Rt (

t)

=

I, det Ri(t)

=

1) describing the rotation of .!fii about .Ai, and ai(t) is the position of .A' ₁ in G( t) (see Fig. 1 ).

()

The velocity of the point.% now reads,

(2.2.2)

Note that since the vectors X and A₁ are vectors in G_0, they do not depend on the time

t.

The vector a₁₍

t)

is the velocity of the point vi _1• The difference X - A₁ can be expressed in the difference x(t) - a₁(t) with the aid of (2.2.1). Then, (2.2.2) can be rewritten as,

(2.2.3)

From the orthogonality of R₁₍ t) it follows that the product

il

1( t) R/( t) is anti-symmetric

(il

₁(t) Ri*(t) = -R₁(t)

ili*(t)).

Denote the product by Ui(t), then there exists with Sl₁₍ t) a vector w₁₍ t) such that for every vector z E IR_3,

(2.2.4)

The vector wi( t) is called the angular velocity of the rigid body .!fi i· With (2.2.4) the equation (2.2.3) can be written as,

(2.2.5).

In order to describe the deformation of the massless deformable body .!fi, the displacement of the point .9'₂ with respect to .9'₁ has to be calculated. Let Pi E G₀ be the position of .9'₁ in the reference configuration. With (2.2.1), the position

Pi( t) E G( t) of S'i at the time treads,

The displacement of the point .9'₂ with respect to the body .i' ₁.can be determined as follows. Consider the mapping x: G(t) 1---+ r(t), defined by,

(2.2.7) x(x,t) = R1*(t)

[x-

~(t)] +Ai'

x

E G(t).

This mapping maps the configuration G(

t)

at the time

t

onto a configuration r(

t)

such that the body .i' ₁ is in its reference configuration. This means that the position vector x( t) of a point$ E .i' ₁ at the time t will be mapped onto the position X of$

in the reference configuration. In particular, for the point .9' ₁ holds,

(2.2.8) P 1 = X(P1( t),t) ·

In general, the bodies .i' ₂ and .i' will not be mapped onto their reference position by

X· This means that the image :r₂₍ t) of p₂₍ t) due to the mapping

x

will differ from P _2. The displacement U( t) of the point .9'₂ with respect to the body .i' ₁ is then given by the difference :r ₂₍ t) - P _2,

(2.2.9)

Using (2.2.8), this equation becomes,

(2.2.10) U(t) =

Ili

*

(t) [p₂(t) - p₁(t)] - (P_{2 -} P₁₎

.

=

Ili

*

(t) y(t) - Y,

where Y

=

P₂ - P₁ and y(t)

=

p₂(t) - p₁(t).

The relative velocity V( t) of the point .9'₂ with respect to the rigid body 2 ₁

equals the time derivative of U( t), since the di placement U( t) is measured from the reference configuration of .i' _1,

(2.2.11)

Multiplying the first term of the right-hand side of (2.2.11) by R₁*(t) R₁(t) (=I), and using (2.2.4), yields,

{2.2.12)

The relative displacement U(t) and velocity V(t), given by (2.2.10) and (2.2.12), respectively, describe the motion of the point !/'₂ of the massless deformable body

~, with respect to the inertial space, if~ is rigidly connected to the inertial space in the point .9'_1.

2.3 The relative orientation a.nd angular velocity of the two rigid bodies

Let.%e~₂ be an arbitrary point of~₂ with position X in G₀ and x(t) in G(t). Its position {( t) in r( t) is given by'

(2.3.1) {(t)

=

x(x(t),t)

=

Ri

*(t)

[x(t) - a₁(t)]

+

A_{1 .}

Using (2.2.1) with i=2, and (2.2.7), (2.3.1) becomes,

(2.3.2) {(t)

=

Ri*(t) Ri(t) (X-A2)

+

x(a.i(t),t) = R(t) (X-A2)

+

ai(t), .

where R(t) = R₁*(t) Ri(t) and a₂(t) = x(~(t),t). The orthogonal mapping R(t)

describes the rotation of the body ~ ₂ from the reference configuration to the configuration r( t), and a₂₍ t) is the position of the point .A₂ in r( t).

{2.3.3)

where 11( t) = R( t) R* (

t)

is an anti-t1ymmetric linear mapping. With 11( t) there exists a vector w( t), the relative angular velocity of.~ ₂ with respect to j9_1, such that

for every vector z E IR_3,

(2.3.4) 11(

t)

z

= w(

t)

>< z .

This relative angular velocity can be expressed in the angular velocities w₁(t) and

w₂₍

t)

of the rigid bodies j9₁ and j9 _2• From the definition of 11( t) and R( t) it follows

that,

(2.3.5)

Together with (2.3.4), this relation implies,

(2.3.6)

hence, the relative angular velocity reads,

(2.3. 7)

In the preceding two sections, the relative displacement, velocity and angular velocity is determined, as well as the relative orientation. The orientation is described by an orthogonal mapping R( t). In the next section, a vector shall be introduced that describes the relative orientation in a linear theory, when the rotations are small. This vector is used in the next two chapters. It shall also be proved that the time derivative of this vector equals the relative angular velocity.

2.4 The relative orientation in a linear theory

The orthogonal rotation tensor R( t), describing the change of orientation of the rigid body ~ ₂ from its reference configuration to the configuration r( t) can be written as, (2.4.1) 00 rpk(t) R(t) =exp (rp(t)N(t)) =

~

- -Nk(t),

£.

k! k=O

where rp( t) is a scalar and N( t) is a real anti-symmetric linear mapping

(N*(t)

= -

N(t)) with llN(t)ll

= 1 (c.f. Kreyszig, 1978, Section

10.5, Theorem 10.5--4, Eqn. 7; rp

=

lliSll and N

=

iS/ rp). Associated with N( t) there exists a vector n( t) such that for every vector x E IR_3,

(2.4.2) N( t) x = n( t) x x.

Then,

II

n( t)

II

= 1 and N( t)n( t)

=

0. The vector n( t) can be interpreted as the axis about which the rigid body ~ ₂ has to be rotated in order to obtain the orientation of ~ ₂ in the configuration r( t) from the orientation of .!fi ₂ in the reference configuration. The number rp(

t)

is the rotation angle ( c.f. Euler's theorem, see Wittenburg, 1977, Section 1). From this it is obvious that the vector n(

t)

is an eigenvector of R( t) with eigenvalue 1. This can also be shown in the following way: using (2.4.2) one can easily derive that the mapping N( t) satisfies,

(2.4.3) N3(t)

+

N(t)

=

0.

Split the series in (2.4.1) into a series containing the odd powers of N(t) and one containing the even powers of N( t), then, with the aid of (2.4.3), one can rewrite (2.4.1) as,

(2.4.4) R(t) =I+ sin ~(t) N(t)

+

(1-cos ~(t)) N2

(t).

From the latter equation it follows that R(

t)n(

t)

=

n(

t).,

hence

n( t)

is indeed an eigenvector of R( t) with eigenvalue 1.

The tensor D(t), defining the angular velocity of.2₂ in r(t), is related to R(t) by D(t)

=

il(t)R*(t). Using the expression (2.4.4) for R(t) leads to,

(2.4.5) ll(

t)

= ~ N +sin ~

N

+ (1 - cos ~) (N

N - N

N) .

When the rotations are small (I~(

t)

I <

1) then (2.4.4) and (2.4.5) reduce to, R(t) =I+ ~(t) N(t),

(2.4.6)

ll( t) = ~( t) N ( t)

+

~( t)

N (

t) .

Now define the vector ~

t),

describing the rotation of .:B ₂ from G ₀ to f ( t), by ~(

t)

= ~( t)n( t). Then, from (2.4.6), it follows that,

R(

t)

x = x

+

~

t)

x x ,

(2.4. 7)

~( t) = ~( t) n( t)

+

~( t)

Ii( t)

=

w( t)

,

in the linear theory.

Now the relative displacement, velocity, orientation and angular velocity is derived, the forces and torques exerting of the rigid bodies can be calculated.

2.5 The forces and torquffi applied to the two rigid bodiffi

At the time

t,

the massless deformable body !iJ is exerted by a force fi( t) in the

point.9'i and a torque mi(t), due to the rigid body!iJi (see Fig. 2).

~(t)

0

Figure 2. The body !iJ exerted by the forces fi(

t)

and torques mi( t).

Since !iJ is massless, these forces and torques satisfy the following two equilibrium

relations,

f1 ( t)

+

f2(

t) = 0 '

(2.5.1)

m₁(t)

+

~(t)

+

y(t) x f₂(t)

=

0.

When the configuration G( t) is mapped onto r(

t)

by the mapping

x,

the loads

(2.5.2)

Multiplying (2.5.1) with R/(t) at the left, yields the following two equilibrium relations for these latter forces and torques,

(2.5.3)

Here, fJ(t)

=

R₁*(t)y(t). The loads Fi(t) and Mi(t) are functions of the relative displacement, the relative velocity, the change of orientation and the relative angular velocity, derived in the preceding sections. Let the force F₂(t) and the torque M_{2 (} t) be given by,

F₂(t)

=

F(U(t),V(t),R(t),w(t))

=

F(t),

(2.5.4)

M₂(t) = M(U(t),V(t),R(t),w(t)) = M(t).

(In the linear theory, the relative orientation is described by the vector ~(t), so in that case the third argument of F and M becomes ~( t), instead of R( t).) The relations (2.5.4) depend on the kind of material the body~ is made of and on the shape of~. They have to be specified for a specific deformable body. From these relations and from (2.5.3), the force F₁(t) and the torque M₁(t) can be expressed in

(2.5.5)

Finally, using the inverse relations of (2.5.2) and the equations (2.5.4) and (2.5.5 ), the forces fi( t) and torques mi(

t),

the rigid bodies exert on the massless deformable body read,

£₁₍ t) = -Ri( t) F(

t) ,

(2.5.6)

m₁(t) = -R₁(t) [M(t)

+

11(t) x F(t)].

~(t)

=

R.i(t) M(t),

The forces

1

i( t) and torques mi( t), which the massless deformable body applies to the rigid bodies, are the opposite of the loads fi(

t)

and mi(

t),

hence,

(2.5.7)

ID1( t) = Ri(

t)

[M(

t)

+

11(

t)

x F( t)] ' ~(t) = -R₁(t) M(t).

2.6 The implementation in MADYMO

In this section, the representations of relative displacement, relative velocity,

change of orientation and relative angular velocity with respect to an inertial

coordinate system E = {E_1, E_2, E_{3 }} are determined. With these representations and

the components of the force F( t) and torque M( t), representations of the forces

f

i(

t)

and torques mi( t) are given with respect to a local coordinate system ei( t) =

{ei(t), aj(t), eMt)}, of the rigid body j5'i, in MADYMO known as the local element coordinate system. These forces are to be added to the equations of motion in MADYMO. The fixed point .Ai of j5' i is chosen coincident with origin of the local

element coordinate system.

el(

t)

. aj(

t)

Let fi( t)

=

{ff( t), fM t), fM t)} be a second coordinate system of the rigid body j3'i, with origin in the point

.91·

The base fi( t) shall be chosen so that in the reference configuration G₀ it coincides with the inertial coordinate system. The base ei(t) shall coincide with a base Ei

=

{Ef, EJ, EH in the reference configuration (see Fig. 3). The latter coordinate system has a fixed orientation to the inertial base E, indicated by the time independent orthogonal tensor S_1,

(2.6.1)

Since Ri( t) represents the rotation of the rigid body ~ _1, the bases fi( t) and E and the bases ei( t) and Ei are related to each other by,

(2.6.2)

Now, MADYMO determines for every rigid body ~i the matrix Ti(t), with respect to the inertial base, of the orthogonal mapping Ti( t) which maps E onto the local element base ei(

t)

of that rigid body,

(2.6.3)

(In MADYMO, T₁(t) is called the cosine direction matrix, see "MADYMO Programmer's Manual", Part II, Chapter 2, Section 2.2.) From the preceding relations it can be derived that the mapping Ri(t) is related to Si and T₁(t) by,

hence, its matrix with respect to the base E reads,

(2.6.5) - 1 R-(t) = - 1 T ·(t)

s ..

-1

Si is the matrix of Si with respect to the inertial base E.

The relative displacement U( t) of the two rigid bodies is given by (2.2.10). If

U( t) is the 3xl matrix made up by the components of U( t) with respect to the

inertial base, and y( t) and Y are these matrices of the vectors y( t) and Y then (2.2.10) yields,

(2.6.6)

u (

t)

=

RI T ( t) y( t) -

y '

In the same way, the representations V( t) and ~(

t)

of the relative velocity V(

t)

and relative angular velocity w( t) with respect to the inertial coordinate ·system can be written as,

(2.6.7)

where. l!i(

t)

and ~i( t) are the representations of the vectors Pi( t) and wi( t),

respectively, with respect to E (c.f. (2.2.12) and (2.3.7)). Since Eis an inertial base the representation of the velocity Pi(

t)

in that base is given by Iii(

t).

The change of orientation of the two rigid bodies is described by the tensor R( t) =

Ri

*( t)Il..i( t). Its matrix R( t) with respect to the inertial base can be expressed in the matrices Ti( t) and Si, by

In a linear theory, the vector ~( t) describes the change of orientation. From equation (2.4.

7)

1 _{it follows that its components ~i( t) with respect to the inertial}

coordinate system read,

(2.6.9)

where Rk₁(t) is the kl-component of the matrix R(

t).

Since the force F ₂₍

t)

and the torque M₂₍ t) are related to U(

t),

V(

t),

R(

t).

and

tiJ(t) by (2.5.4), their representations

E

₂(t) and M₂(t) with respect to the inertial bruie are related to U ( t), V (

t),

R( t) and ~(

t)

by,

F2(t) = F(V(t),V(t),R(t),~(t)) = F(t), (2.6.10)

M₂(t) = M(V(t),V(t),R(t),~(t)) = M(t),

(in the linear theory the matrix

R(t)

in (2.6.10) is replaced by the matrix ~(t)

=

[~

₁

(t), ~

₂

(t), ~

₃

(t)]T). Using the relations (2.6.3), (2.6.4) and (2.5.7), the representations

L(t)

and mi(t) of the loads fi(t) and Iili(t), respectively, with respect to the local element coordinate system ei(

t)

can be expressed in F(

t)

and

M(t) by,

m

1( t)

=

s.

1

[M(

t)

+

11(

t) x

F(

t)] ,

(2.6.10)

These latter quantities have to be added to the equations of motion in MADYMO.

2. 7 Input description

The input data must be placed at the end of the input file "DATA.DAT" after

the keyword USER INPUT, but before the keyword END INPUT. The input to be

given is described below.·

USER INPUT

'MASSLESS DEFORMABLE BODY'

SYS IPAR EL PARl

x

y

PAR2 Z TEXT PAR9

(These two data records have to be repeated twice.)

YX yy

YZ

Variable Type Units Definition

SYS Integer System number of rigid body .!fi i to which the

Variable Type Units Definition EL

x

y

z

IPAR PARl PAR2 PAR9 TEXT YX

yy

YZ

Integer Real m Integer Real Character -Real m

Element number of rigid body$ i to which the massless deformable body in connected

The coordinates of the point .9'i in which the massless deformable body is connected to$ i. The coordinates are expressed in the local coordinate system if$ i·

Parameter for orientation definition of the

coordinate system in the point .9'i with respect to the local coordinate system of$ i:

IP AR= 1: successive rotations IP AR = 2: screw axis method IP AR = 3: direction cosine

(c.f. "MADYMO User's Manual", Chapter 4, Section 4.8)

orientation specification

Identifier of the element EL of system SYS (max. char. = 15)

The coordinates of the point .9'₂ with respect to .9'₁ in the reference configuration. The coordinates are expressed in the inertial coordinate system

After this, the deformable body characteristics have to be specified. These characteristics depend on the kind of material. In the next two chapters a linear elastic beam and a linear visco--€lastic beam are modelled. The required input data for these beams are given in these chapters.

3

A

linear elastic beam with a circular cross-section

3.1 Introduction

As an example of how to use the subroutine DEBOFO, a linear elastic beam with a circular cross-section shall be implemented in MADYMO. This beam is rigidly connected to two rigid bodies in the centres of its final cross-sections. In order to use DEBOFO, the force and torque exerting on one final cross-section, while the other is rigidly connected to the inertial space, have to be derived as functions of the displacement, the velocity, the change of orientation and the angular velocity of that final cross-section. Since the beam is made of an elastic material, these loads are independent of the velocity and the angular velocity.

3.2 The loads on the linear elastic beam

Consider a linear elastic beam, length

£

_{0 ,} with a circular cross-section, radius R.

This beam is rigidly connected to the inertial space in one final cross-section (&' _1). On the other end cross-section (&' ₂₎ a torque M and a force F through its centre are applied (see Fig. 4). Let E = {E_1, E_2, E₃} be an inertial coordinate system with

origin in the centre of&' _1. The E₂ axis of this base coincides with the axis of the beam in the undeformed configuration, which here is taken as the reference configuration. In the preceding chapter the diplacement U = UkEk of the centre of

&' ₂ has been determined, as well as the tensor R, characterizing the change of orientation of the cross-section &' _2• Since the beam in made of a linear elastic material, all deformations are assumed to be small. Thus, the vector 'P can be used to describe the change of orientation. Its components 'Pk with respect to E can then

f,R,E,G F(t)

~M(t)

Figure 4. A linear elastic beam exerted by a force F(

t)

and a torque M(

t).

be seen as the angles about which the cross-section rfl' ₂ rotates: ~1 and ~3 are the

bending angles about the E₁ and E₃ axis, respectively, and ~_{2 is the torsion angle of}

the beam.

Let F

=

FkEk and M

=

MkEk be the representations of F and M, respectively,

with respect to E. From the Theory of Elasticity it is known that the extension U₂ is related to the normal force F₂ by,

(3.2.1)

where Eis Young's modulus and S = -xR2 is the area of the cross-section of the beam. The angle ~_{2 can be expressed in the torque M2} causing the torsion by,

(3.2.2)

Here,

G

is the shear modulus and IP the polar moment of inertia of the cross-section, defined by,

(3.2.3)

IP =

f

_(x12

+

X:J2) d<T =

~

R4 .

s

force F_3• The relation between the displacement U₃ and the angle y?₁ and these loads is given by, (3.2.4) Mi fo2 F3fo3 U3= + _ _ 2EI 3EI Milo F3lo2 Y'1 = - - + _ _ EI 2EI

Here, I is the moment of inertia of the cross-section about the E₁ axis or the E₃ axis, defined by,

(3.2.5) I=

f

~2

dq

=

f

X12 dq

=

f

R4.

s

s

Analogous equations hold for the bending of the beam about the E₃ axis,

(3.2.6) Mafo2 F1foa Ui=-

+

-2EI 3EI M3lo F1 fo2 y ? 3 = -EI 2EI

Inverting the equations (3.2.1), (3.2.2), (3.2.4) and (3.2.6), the forces Fk and torques Mk can now be expressed in the displacements Uk and the angles Y'k· This results in,

(3.2. 7) 12£/ 1 Fi =

(3 (

U1

+

2

f o

~3)

'

0 12£/ 1 F3

=

(3 (

U3 -

2

f o

~1)

'

0 12E/ ( l l ) Mi =

(3 -

2 U3

+

3 f o ~t , 0

With these relations the subroutine LEBEFO can calculate the forces and torques that cause the linear elastic beam to deform. These forces are passed to the subroutine DEBOFO, which adds them to the equations of motion.

The elastic beam can be connected to two MADYMO elements by placing the keyword 'ELASTIC BEAM' after the vector Y has been specified in the MADYMO input file. Next, the radius R (in m), Young's modulus E and the shear modulus G (both in N /m2) have to be given.

3.3 Some simulations

been performed, that can easily be checked in an analytical way. Both stretching, torsion and bending of the beam were simulated. In the stretching and torsion simulations, the beam was connected to the centres of two identical spheres. In the case of stretching, these centres were initially placed at a distance of £₀

+

8,

without initial velocities. The length of the beam is then a harmonic function of the time and the period of the vibration reads,

(3.3.1) _{T= 2:r (} mf₂

_~S

0

)~

_'

where mis the mass of one sphere. The quantities have been chosen such that the period equals 1 s. The simulation performed by MADYMO resulted in the graph presented in Fig. 5. This graph gives the length of the beam as a function of the time. The dotted line indicates the underformed length of the beam which is the equilibrium position of the vibration. One can see that the period of vibration indeed equals 1 s. ,... E ...,, QJ u ffi 1 .005 ,..., .\!! 0 1. 000 0.995 0.990-+-~----.-~~-.-~----.-~~-.-~--. 0 2000 4000 Time(ms)

In the case of the torsion simulation, the two spheres were initially twisted around the axis of the beam about an angle

o/2,

such that the two final cross-section of the beam were twisted about an angle

o.

Again, no initial velocities and angular velocities were given to the spheres. The relative twist of the two spheres is then a harmonic function of the time, with a period of vibration, equal to,

(3.3.2) ₍ mr

2

l0

)!

T = 21' .

5G/P

Here, r is the radius of one sphere. As in the preceding simulation, the quantities were chosen such that the period equaled 1 s. After the simulation that MADYMO performed, the sine of the relative twist of the two spheres was drawn as a function of the time, resulting in the graph of Fig. 6. Again, the dotted line indicates the equilibrium position of the twist, which is 0 rad. From Fig. 6 one can see that the period of the vibration indeed equals 1 s.

oJ c QJ E QJ u

m

a 0.04 0.02 .'.!! -0.00 u a

5

-0. 02 u I

x

-0.04 0 2000 4000 Time( ms)

In the bending simulations one end cross-section of the beam was connected to the centre of a very heavy sphere while the other was connected to a light sphere. The inertia of the heavy sphere were taken so large that its motion can be neglected. Then, one can assume that the beam is connected to the inertial space in one final cross-section. The centres of the two spheres were initially placed on the E₂ axis of the inertial coordinate system and the light sphere was given an initial velocity of v₀E₃ and angular velocity w₀E_1. The component U₃ of the centre of the light sphere along the E₃ axis is then given by,

(3.3.3)

Here, Ai and Bi are constants depending upon the initial conditions. The constants

wi depend upon the properties of the beam (E, I, l_{0 )} and upon the mass m and the

radius r of the light sphere. An analogous equation holds for the rotation angle ~_{1 of}

the light sphere about the E₁ axis. From (3.3.3) it can be seen that U₃ is a superposition of two vibrations. The periods of these vibrations read Ti = 2;r / "'i·

In the simulation, the initial velocity and angular velocity were chosen such that only one vibration occured (i.e. A₂ = B₂ = 0). The properties of the beam and of the light sphere were determined such that the period of the remaining vibration equaled 1 s. In Fig. 7, the component U₃ is drawn as a function of the -time. The dotted line again indicates the equilibrium position of the component (i.e. U₃ = 0). From Fig. 7, it can be seen that the period of the vibration indeed equals 1 s.

,_J c QI E Q/ u lil ~ Q. 0.002 0.001 .~ -0. 000 D Q. E 0

y

-0 .001 N -0.002-+-~--i.--~-,.-~~...-~--,.-~----, 0 2000 4000 Tlme(rns)

4 A linear visco-€lastic beam with circular cross-section

4.1 Introduction. Basic equations

As a second example of how to use the subroutine DEBOFO, a model shall be made of a massless linear viscC>-€lastic beam with a circular cross-section. For this, the relation between the forces and torques on one hand and the displacements and rotations on the other, has to be derived. The deformations are assumed to be small, such that a linear theory can be used. In that case, the rotation angles, as derived in Chapter 2, Section 2.4, are the bending and torsion angles.

Consider a linear visco-elastic beam, length £_0, with a circular cross-section, radius R. Let one end of the beam be rigidly connected to the inertial space, while the other is exerted by a force F( t)

=

Fi( t)Ei and a torque M( t)

=

Mi( t)Ei (see Fig. 8).

~~E-3_•_E_2~~~f-,R~~~~~~)

M(t)

Figure 8. A visco-elastic beam exerted by a force F(

t)

and a torque

M(

t).

The unknown quantities in this problem are the displacements ui(x,t) and the stresses t₁j(:x,t) of a point.% of the beam with position x ( ui(x,t)

=

(u,Ei), tij(x,t)

=

(Ei,TEj), where u(x,t) is the displacement vector and T(x,t) the stress-tensor). Define as auxilary quantities the linearized deformations eij(x,t) and the deviatoric

deformations € ij(x,t) and deviatoric stresses si;(x,t),

(4.1.1)

where 'i denotes the partial derivative

8/ 8xj

(x

=

xjEj) and b'ij is the Kronecker-delta.

Since the beam is massless the stresses ti;(x,t) have to satisfy the following equilibrium equations,

(4.1.2)

In addition to these three equations, six equations that describe the behaviour of the material are needed to determine the solution. (Note that there are nine unknown quantities: three displacements and six stresses, since the stress-tensor is symmetric.)

4.2 Constitutive equations

Usually, compression of a visc<>-€lastic material is assumed to be elastic (in general the changes in shape are much larger than the changes in volume). This yields the following relation bet ween the trace of the linearized deformation tensor and the trace of the stress tensor,

where I< is the compression modulus.

The relation between fij(x,t) and silx,t) (shear) is given by the general constitutive equation for a linear visco-elastic material,

(4.2.2)

The integral on the right-hand side is called the hereditary integral. The function

~ G ( t) is the creep function for shear. This function represents the respons of the material on a load, which is 0 for t

<

0 and 1 for t

>

0. In Section 4.6 a specific function shall be chosen for ~G(t): the creep function for a standard-linear solid. First the boundary conditions have to be specified.

4.3 Boundary conditions

To the equations ( 4.1.l ), ( 4.1.2) and ( 4.2.2) boundary conditions have to be added (six in each point of the boundary). Since the beam is rigidly connected to the inertial space in

x.i

=

0, the displacement has to be zero there,

(4.3.1)

Next, the part of the surface x₁2

+

x₃2 = R2 of the beam has to be stress-free. This results in three relations on de edge of the circle. Two of these are trivially satisfied, while the third reads,

(4.3.2) _{~1 n1}

₊

_{~2 ~}

=

_{0 ,}

On the end cross-section Xi = l₀ a force F( t) and a torque M(

t)

is exterted, which yields, F( t)

=

f

ti

₂E₁ dq , Xi

=

l_{0 ,}

s

( 4.3.3) M(t) =

J

(xx

ti

₂Ei) dq, Xi= l_{0 •}

s

With these boundary conditions the solution is unique.

4.4 The transformation of the problem

The problem stated in the preceding sections can be solved by transforming it by a Laplace transform with respect to the time

t.

Suppose every quantity A( t) is defined for

t >

0, and let A( t) = 0 for t

< 0.

Define the Laplace transform A( s) of

A(t),

by

(4.4.1)

00

A(s)

=

.t'{A(t)}

=

J

A(t) e-st dt.

0

The Laplace transform has, among others, the following property,

( 4.4.2)

.t'{A(t)}

=

s.t'{A(t)} - A(O)

=

s A(s) - A(O) ,

(c.f. Kreyszig, 1988, Chapter 5, Section 5.2, Theorem 1). Next, for the product of two Laplace transforms the convolution theorem holds,

( 4.4.3)

t

..t{f A(t- r) B(r) dr} =..t{A(t)} ..t{B(t)}

=

A(s) B(s),

0

( c.f. Kreyszig, 1988, Chapter 5, Section 5.6, Theorem 1 ).

Due to the Laplace transform, the equations (4.1.1), (4.1.2) and (4.2.1) and the boundary conditions ( 4.3.1), ( 4.3.2) and ( 4.3.3) are transformed into analogous equations for ui(x,s), eij(x,s), tij(x,s), sij(x,s), tij(x,s), .Fi(s) and Mi(s). Using the properties ( 4.4.2) and ( 4.4.3), equation ( 4.2.2) is transformed into,

( 4.4.4) with, ( 4.4.5) 1 G(s)

=-_-s?G ( s)

Now the problem is transformed, the solution can be obtained easily.

4.5 The solution of the problem

Note that the problem made up by the transformed equations ( 4.1.1), ( 4.1.2), (4.2.1), (4.3.1) - (4.3.3) and equation (4.4.4), fits with the static linear elastic problem solved in Chapter 3, where /(and G(

s)

are the compression and the shear modulus, respectively. The transformed forces and torques applied to the end cross-section of the beam can be expressed in the displacements Ui(s) = ui(0,(₀,0,s)

of the centre of the end cross-section and the rotation-angles ~i( s) of that cross-section, by

_ 12.E(s)/ [- ₁ _

J

F1(s) = l ₃ U1(s)

+

₂ l0 ~

3

(s)

,

0 _ 12.E(s)I[_ ₁ _ ] F3(s) = l 3 U3(s) -

2

£0 ~1(s) , 0 ( 4.5.1) _ 12E( s)I [ l _ l _ ] M1(s)

=

-2

U3(s)

+

3

l0 ~1(s) , f, 2 0 _ _{12.E(s)/ [ 1 _ 1} _

J

.

M₃(s) = l 2 2 U1(s)

+

3 l0 ~

3

(s) , 0

Here, I is the moment of inertia of the cross-section S, IP is the polar moment of inertia of S, and Sis its area,

-(4.5.2) I =

f

~

2 d<T =

f

Xi 2 d<T

=

i

R4 '

s

s

I -

f

(x 2

₊

_x 2 ) d<T - ! R4 p - 1 3 -2 '

s

S =

f

d<T

=

r R2 ,

s

and E( s) is Young's modulus, related to [(and G(s) by,

(4.5.3)

9KG(s)

E(s) = _

3K

+

G(s)

The solution of the original problem can now be found by determining the inverse Laplace transform of the relation ( 4.5.1 ), using the convolution theorem (4.4.3). Before doing this, a choice shall be made for the creep function ,G(t). The solution shall be derived for a standard linear solid. In the next chapter, this and two other simple visco--elastic models shall be discussed.

4.6 Some simple visco--elastic models

The simplest visco--elastic models, the Maxwell model and the Kelvin-Voigt model, can be presented as a series connection and a parallel connection, respectively, of a linear spring and a linear viscous damper (see Fig. 9).

T/ c F " - - - + - - - 1 11(

t)

- - + €( t)

(a)

1 - - - - t 11( t) - - +

f( t)

(b)

The creep function ~( t) of such a model can be found by exerting it by a force 11( t) = H( t), the Heaviside function. The function ~( t) then equals the extension

€( t). A more complicated model is the standard linear model. This model can be represented by a series connection of a Kelvin-Voigt model and a linear spring (see Fig. 10).

fL----:?i+---+ 11( t)

----+ €( t)

Figure 10. The standard linear model.

Note that for c₁=0 the Maxwell model arises from the standard linear model, and that the Kelvin-Voigt model is obtained by taking the limit for ~to infinity.

Now, the creep function of this standard linear solid shall be derived. Let _{f 1}(t)

be the extension of point ~ (i.e. the extension of the spring .A.!iJ and the displacement in the damper ~!!J) and let €₂(t) be the extension of the spring ~.5f.

The total displacement €( t) of .5f, due to the force D"( t) then equals their sum,

(4.6.1)

The force in the spring~ .5f and the total force in the parallel connection .A.:IM !lJ are both equal to O'( t), so for the spring ~.5f holds,

(4.6.2)

force in the spring .A~ (u_1(t) = c₁t_1(t)) and u_2(t) is the force in the damper'& 9J

(u₂(t) = 77€₁(t)), hence,

Using the equations (4.6.1) and (4.6.2), _{t 1}(t) can be expressed in u(t) and t(t) by,

( 4.6.4)

Substituting this relation into ( 4.6.3) finally yields,

( 4.6.5)

with /i

=

cJ17.

Solving t:( t) from the linear ordinary differential equation of first order leads to,

t

-

t

[

1 fl

+

/2 . ] T

(4.6.6) t:(t)=e/1 J.-u(r)+ u(r) eli dr.

-00 77/2 77 / 2

By definition, the creep-function ¢(

t)

equals the extension

t:( t)

if u(

t)

= H(

t).

Then,

'1( t) =

c5(

t), the Dirac-delta function, and,

t

t [ 1 /1

+

/2 ]

(4.6.7) ¢(t) = e-ft J - b(r)

+

H(r) eliT dr.

-00 77/2 77 / 2

( 4.6.8)

while for t

<

0, ;( t)

=

O and

;(o)

is defined by,

(4.6.9) _{;(o) :=}

u~

_{;(t) =}

7/~2.

As stated before, the Maxwell model arises from the standard linear model by taking c₁ = 'l/1i

=

0. Using the Taylor expansion of the exponential function, from ( 4.6.8) it follows that,

(4.6.10)

hence the creep function of a Maxwell model is given by,

( 4.6.11) 1 t ;( t) = 1112

+ 7;

I t

>

0 ! 1

;(o)

= 1/12 •

Similarly, one can obtain the Kelvin-Voigt model from the standard linear solid, by taking the limit

c.i

= 111₂ ~ oo. Then, ( 4.6.8) and ( 4.6.9) yield,

( 4.6.12)

In the next sections, the solution, as derived in Section 4.5, shall be worked out further for a standard linear model. Hence, the creep function for shear, ~ G ( t), shall be taken equal to expression ( 4.6.8).

4. 7 The solution for the standard linear solid

In order to determine the solution, the Laplace transform ~G(s) of the creep function is needed. Using the definition (4.4.1), one obtains,

(4.7.1)

00

1

+

12

+

s ~ (s)

=

J

(-1

+

_1 (1 - e-11t)) e-st dt

=

_ 1 _ _ _ _

G 0 . 1/"f2 1/i1 1/i2 s (s

+

ii)

The shear modulus and Young's modulus are then given by,

( 4.7.2) with, (4.7.3) _ 1 [ G0

l

G(s) = _ = G₀ 1 - , s~G ( s) 1111

+

G0

+

17s _ 9KG(

s)

[

!Eo

l

E(s) = _ = E₀ 1 -3[(

+

G( s) 1/ii

+

_!E0

+

1JS ' 9KG₀

Eo=----3[(

+

G₀

The solution is obtained by determining the inverse Laplace transform of the equations (4.5.1). These equations all have the same form,

(4.7.4) L( - s)

=

C₀ [ 1 - 1 -

11] -

d( s) , 1

+

s

where L represents the force or the torque, C₀ is an elasticity constant ( C₀ = G₀ or C₀

=

E_{0 )} and 1

=

1₁

+

G₀/71 or 1

=

1₁

+

E₀/(371). The function d is a linear combination of a displacement ui and an angle ~j· The inverse Laplace transform of the latter equation is easy to determine with the aid of the convolution theorem (4.4.3). The inverse Laplace transform of (1+st1 _{reads e-1t, hence,·}

(4.7.5)

t

=Co d(t)- C

0 (1-11)

J

d(r) e-1(t- r) dr.

0

In ( 4. 7.5), L( t) is expressed in the deformation d( t) under the condition that

L(t) = O, d(t) = O, fort~ 0 (i.e. Fi(t)

=

0, Mi(t) = O, ui(t)

=

0, ~i(t) = 0, fort~ 0). One can also determine L(

t)

if at time

t

₀

>

0 the load and the deformation are known. For this purpose, rewrite (4.7.5) as

(4.7.6)

t

₀

t

L(t)

=

C₀ d(t)- C₀ (1-1₁)

(J

+

J)

d(r) e-1(t- r) dr.

0

t

₀

Using (4.7.5), the integral from 0 to

t

₀ can be expressed in L₀ = L(t_{0 )} and d₀

= d(t

_0). Substituting the obtained expression into (4.7.6) yields,

(4.7.7)

t

- Co ( r- r

1) e-r(t- to)

f

d(r) e-r(r - to) dr.

to

Finally, the latter equation can be rewritten as,

(4.7.8)

t

- C

0 ( r - r 1) e -r( t - to)

J (

d( r) - d0) e -r( T - to) dr .

to

Note that the left-hand side the increase of the load represents in the time interval

[ t, t_0] and that it is expressed in that time interval and the increase of the

deformation d( t) - d₀ in [ t, t_0].

4.8 The formulation in ODE's

The equations ( 4. 7. 7) and ( 4. 7.8) are integral representations of the load

L(

t).

For calculations it is often more convenient to have the relation between the load and the deformation expressed as an ordinary differential equation (ODE).

Differentiating ( 4. 7. 7) with respect to the time

t,

yields

(4.8.1)

t

+Co r (r- r₁) e-r(t- to)

J

d(r) e-r(r - to) dr

The integral on the right-hand side of ( 4.8.1) can now be expressed in L(

t),

d(

t)

by using (4.7.7). Substituting the result in (4.8.1) finally results in the following ODE for Land d,

(4.~.2)

L

+

r

L = C_{0 (}

d

+

r

₁ d) .

Together with the initial conditions L(

t

_{0 )} = L_{0 ,} d(

t

_{0 )} = d₀ the solution of ( 4.8.2) exists and is unique.

Writing out ( 4.8.2) for the various equations ( 4.5.1 ), one obtains the following six ODE's describing the relation between the forces and torques on one hand and the displacements and rotation angles on the other,

with initial conditions,

( 4.8.4)

Here, /E

=

11

+

E0/(3TJ) and /G

=

11

+

G0/TJ.

As stated before, the Maxwell model can be derived from the standard linear

solid by taking c₁ = 0, hence, for a Maxwell material the equations ( 4.8.3) hold

with 1₁

=

0. The Kelvin-Voigt material arises if c_{2 -+} oo. For E₀ this yields,

(4.8.5)

9Kc₂

E₀

=

-+ 9 [( , c_{2 -+} oo ,

3[(

+

c₂

while G₀ = c₂-+ oo. This means that problems arise in equation ( 4.8.3)5• However,

deviding this equation by G₀ gives,

(4.8.6)

Now let G₀ -+ oo, then from ( 4.8.6) it follows that,

(4.8.7)

This means that for a Kelvin-Voigt material, the forces and torques are related to

the displacements and rotation angles by the equations ( 4.8.3)1 - ( 4.8.3)4, ( 4.8.3)6

4.9 The implementation in MADYMO

For the implementation in MADYMO, the beam is connected in its end cross-sections to two rigid bodies. The subroutine DEBOFO determines the relative displacements and rotation angles of these two rigid bodies. Then the ODE's ( 4.8.3) describe the forces and torques that have to be passed to DEBOFO. These differential equations are integrated during the simulation by a subroutine VEBEFO. This is done by using the same methods as MADYMO: a 4th order Runge-Kutta integration method with a fixed time step, or a 5th order Runge-Kutta Merson method with variable time step.

The input required for the visco-elastic beam has to be placed in the input file after the vector Y is specified. First the keyword 'VISCO-ELASTIC BEAM' has to be given. After that an integer MEDIUM, indicating the medium the beam is made of, must follow (MEDIUM

=

1: Maxwell-material, MEDIUM

=

2: Kelvin-Voigt material, MEDIUM

=

3: Standard linear solid). After that, the radius R (in m), the

compression modulus K (in N/m2_{), the damping coefficient 1J (in Ns/m'l.) and the} two spring constants c₁ and c₂ (in N /m2) must be given. If the medium requires only one spring constant (Maxwell, or Kelvin-Voigt) Ci is not used. Finally the

pre-tensions FiO (in N) and MiO (in Nm) must be given.

4.10 Some simulations

In order to test the subroutine VEBEFO, some simulations were performed that can be checked in an analytical way. Both stretching and torsion of the beam were simulated. As for the linear elastic beam, the visco-elastic beam was connected in its end cross-sections to the centres of two indentical spheres. Both the centres of the spheres were initially placed at the E₂ axis of the inertial coordinate system at a distance

£

₀ of each other, such that the beam was in its undeformed configuration (the reference configuration), and no initial forces or torques were exerted on the

beam. In the case of the stretching of the beam, the spheres were given initial

velocities along the E₂ axis that were equal in magnitude but opposite in direction.

Applying Newton's Law to the two spheres, yields that the force F( t)

=

F_{2 (} t)E₂ that stretches the beam, is related to the length l( t) of the beam by,

(4.10.1) F₂(t) = - -1 <Pl

2

m - .

dt2

Equation ( 4.8.3)2 describes the relation between the force F₂₍ t) and the extension

U₂₍ t)

=

l( t) - l₀ of the beam. Substituting ( 4.10.1) into this equation yields the following 0 DE for the length

l (

t) of the rod at the time t,

<f3l <Pl dl

(4.10.2)

_-

₊

_{IE -}

₊

a2 -

+ /

1 a2

l

= / 1 a2

lo '

dt3 _dt2 _dt

where a2 = 2E₀S/(ml_{0 ).} A particular solution of this ODE is the constant solution

l(t) = l_0• Substituting l(t)

=

e,\t in the homegeneous equation gives the characteristic equation of ( 4.10.2),

(4.10.3) ,\ 3

+

I ,\ 2

+

a2 ,\

+ /

1 a2 = 0 .

E

-The roots of this cubic equation are not easy to find in an analytical way. However,

for a Maxwell material ( 1₁

=

0, see Sec. 4.8), ( 4.10.3) reduces to,

( 4.10.4) ,\ ( ,\ 2

+

IE ,\

+

a2) = 0 '

hence, ,\₀ = 0 is one solution. The other two roots depend on the sign of D =

(4.10.5)

~ _{1, 2} = -; E/2 ± iw , D = - w2

<

0 ,

The latter relations imply that the general solution of ( 4.10.2) reads,

f( t) = A + exp (-;E/2) ( B cosh (/)t)

+

C sinh

(/Jt)) ,

D

>

0 ,

(4.10.5) f(t)

=A+

exp (-;E/2)

(n

+Ct),

D=

0,

f(t)=A+exp(-rE/2)

(ncos(wt)+

Csin(wt)),

D<O.

Here, A, B and C are constants, depending on the initial conditions. In the case

D >

0 is system is overdamped, if

D

<

0 the system is underdamped and if

D

=

0 the system is critical damped. Note that if the system is underdamped, the solution is a damped vibration with period T = 2;r /

w.

Both overdamping and underdamping were simulated for a Maxwell medium. In the case of an underdamped system, the material properties were chosen such that the period of the damped vibration equaled 1 s. In Fig. 11, the calculated length of the beam is drawn as a function of the time in the case of an underdamped system and an overdamped system. The dotted lines in these graphs indicate the stationary solution, lti m f(

t)

= A. From Fig. 11 it can be seen that the period of the vibration

'-+oo indeed equals 1 s.

EnGTH ,... 2.30 E OJ _2.25 u c rn ,.., V1 _2.20 0 2. 15 2. 10 2.05 2.00 D 2000 (a) 4000 Time(ms) LEnGTH 3.2

---3.0 2.8 2.6 2 . .:! 2.2 2.0-t-~--.-~--,r--~-.--~--.-~---. 0 2000 {b) 4000 Time(ms)

Figure 11. The length of the beam. Underdamped (a) and overdamped {b).

In the case of torsion of the beam, the spheres were given an initial angular

velocity along the

Bi

axis, equal in magnitude, but opposite in direction. In the same way as above, now using the Law of Moment of Momentum and equation

( 4.8.3)5 one can derive the following 0 DE for the relative twist ~

2

( t) of the two spheres,

(4.10.5)

where a2 = 5G₀

Ip/(mr2f

_0). This is an analogous ODE as (4.10.2) and for a Maxwell model the general solutions are of the same form as ( 4.10.5). In this case only the underdamped system was simulated, were the material properties were chosen such

that the period of the damped vibration equaled 1 s. In Fig. 12 the sine of the

relative twist of the two spheres is drawn. The dotted line indicates the stationary

,... 1 . o E ' - '

...

c _0.8 ru E ru u m _0.6 a. t.ll 0 a. 0.4 E 0 u I 0.2

x

0.0 0 2000 4000 Time( ms)

Figure 12. The sine of the relative twist of the two spheres for a Maxwell model.

In case of a Kelvin-Voigt model, instead of using (4.8.3)5 one has to use the equation ( 4.8. 7) for the relation bet ween the relative twist and the torque. Then a second order ODE arises that can be solved in an analytical way. _Again three types of solutions are obtained: underdamping, critical damping and overdamping. Only the underdamped system was simulated such that the period of the damped vibration equaled 1 s. The sine of the relative twist of the two spheres is drawn in Fig. 13. The dotted line again indicates the stationary solution ~₂(

t)

= 0. From· this figure it can be seen that the period again equals 1 s.

,... 1 . 0 E .._,, ..., c ru E OJ 0.5 u m 0. Ul D 0. E 0 u I

x

-0.S 0 2000 4000 Tirne(rns)

5 Conclusions and suggestions for further improvements

In Chapter 2, the user-defined subroutine DEBOFO is developed. This subroutine brings into account the forces and torques a massless deformable exerts on the two rigid bodies to which it is connected. In order to use this routine a user has to write a subroutine in which the forces and torques that apply to one point of connection of the massless deformable body, while the other is connected to the inertial space, are calculated. In the Chapters 3 and 4 this is done for a linear elastic beam (the subroutine LEBEFO) and for a linear visco-Blastic beam (the subroutine VEBEFO), respectively.

The subroutine DEBOFO can now handle only one massless deformable body that is connected to two rigid bodies. When the subroutine is implemented in MADYMO, the possibility of connecting any two rigid bodies by a massless deformable body is needed. This can be done by splitting DEBOFO into two parts; one part that calculates the relative displacement, velocity, change of orientation and angular velocity of two rigid bodies and one part that calculates the forces and torques on these rigid bodies. A user can then call these routines in his own user-defined routines that calculate the forces and torques that exert on the massless deformable body he wants to implement.

The subroutine can be improved further by adding the possibilty of connecting a massless deformable body ii:J. one point to the inertial space, or to a null system. This requires only minor changes in DEBOFO, since the program copies the data needed to determine the relative motion of the two rigid bodies, into local array's and it uses the MADYMO library routines TRANS and VELACC (see "MADYMO Programmer's Manual", Part II, Chapter 4).

References

Kreyszig, E., 1978, "Introductory Functional Analysis with Applications", John Wiley & Sons. Inc., New York.

Kreyszig, E., 1988, "Advanced Engineering Mathematics", 6th edition, John Wiley

& Sons. Inc., New York.

"MADYMO Programmer's Manual", October 1988, Report TNO Road Vehicles Research Institute, Delft.

"MADYMO User's Manual 3D", October 1988, Report TNO Road Vehicles Research Institute, Delft.

Vosbeek, P.H.J., "Two Contributions to the Modelling of the Human Body in MADYMO", August 1989, Master's thesis, Eindhoven University of Technology, Eindhoven.

Wittenburg, J., 1977, "Dynamics of Systems of Rigid Bodies", B. G. Teubner, Stuttgart.