A detailed elaboration of problem 1 : a rigid sphere in a rigid cylinder, having point contact and interconnected by nonlinear elastic line elements

A detailed elaboration of problem 1 : a rigid sphere in a rigid

cylinder, having point contact and interconnected by nonlinear

elastic line elements

Citation for published version (APA):

Sauren, A. A. H. J. (1985). A detailed elaboration of problem 1 : a rigid sphere in a rigid cylinder, having point contact and interconnected by nonlinear elastic line elements. (DCT rapporten; Vol. 1985.051). Technische Hogeschool Eindhoven.

Document status and date: Published: 01/01/1985

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

providing details and we will investigate your claim.

TO W E 83-05 ADDENDUM

A DETAILED ELABORATION OF PROBLEM 1 :

a rigid sphere in a rigid cylinder, having point contact and interconnected

by nonlinear elastic line elements.

WFW 85.05 1

Contents

Introduction 3

PART I

AHALYTICAL TREATMENT O F THE PROBLEM 1.1

1.2 The rotation tensor

I.3 Formulation of the problem 1.4

Geometrical description o f the cylinder and the sphere

Analytical solution to the problem 1 . 4 . 1 Contact conditions

1.4.2 The external loads on the sphere

1.4.3 The loads on the sphere due to point contact

1.4.4 The loads on the sphere due to kinematic constraints 1.4.5 Loads on the sphere due to the line elements

1.4.6 Force and moment equilibrium

PART I1

EJUMERICAL T ~ E A T ~ ~ ~ T OF THE PROBEER

11.1

Determination of ötp, O$ and öw

11.2 Listings of user written subroutines 11.3 Input data 3 7 9 10 10 12 12 12 13 15 17 20 3 6

Introduction

In this report an elaborate version is presented o f example 1, given in report W E 83-05. The analytical and numerical treatment of the problem will be discussed in Part I and 11, respectively. The reader is supposed to be familiar with the theory presented in report WE 83-05. The numerical

solution of the problem is obtained with the software that is described in WE 83-05. Double numbered equations in the following refer to equations in WE 83-05.

The configuration o f the problem is constituted by a rigid sphere (radius b) within a rigid hollow cylinder (radius a). There is point contact between the two bodies that are connected by three identical line elements (Fig.1).

x

For reasons of simplicity the problem is formulated in such a way that there will be only motions in the x-z-plane, i.e. the problem is two-dimensional. In the analytical treatment this will be used whereas the numerical analysis will be applied to the full three-dimensional problem.

PART I

ANALYTICAL TREATMENT OF THE PROBLEM

1.1

Geometrical description of the cylinder and the sphere.

For the description o f the cylinder we will use the Cartesian xya-coordinate system and the orthonormal vector base e -, =[ex e e

IT

as shown in Fig.2. The

polar coordinates q~ and

z

will serve as surface parameters (see Fig.2). The position vector c of a point on the cylinder is written in terms o f its

-9 - @ + +

Y Z

-@

components with respect to the vector base as

.r

d 4 4 4

c = a cos(p)e t a sin(p)e t z eZ.

X Y

For the tangent vectors to the parametric curves ( z = constant and rp =

constant) holds :

d 4

c _z = e z .

By use of the definition

4 d

(31

the unit normal vector on the cylindrical surface is written as

( 5 ) d d d n =

-

cos(rp)ex - sinfq)ey ~ ~ ~~ ~~~~~

where the choice sn =

-1

has been made so that equation ( 5 ) represents the unit inward normal vector on the cylinder.

Describing the geometry o f the sphere we will use the Cartesian ap-y-

coordinate system, having its origin O Pn the center o f the sphere, and the A

Fig. 3

+ + " b

orthonormal vector base E = E

:

IT

as depicted in Fig. 3 . The spherical

c

$

r

coordinates cp and 8 (see Fig. 3 ) will serve as surface parameters. In the reference configuration of the system holds, by definition, for the vector

- e * +

ba-sec e = E = E and for ~~ the ~ orthonormal -~ rotation tensor R = Ro= - I, ~~ yielding

~ ~ - ~ - ~

c c ,o

for the position vector

c

of a point on the sphere in the reference configuration :

-3 9

= b cos(8)cos(cplex

+

b cos(B)sin(p)e t b sin (@)ez. Y

A

It is noted that the components of

c

with respect to

:

are constants.

Changes of the orientation of the sphere (and thus of

6)

A, a.re described by R.

The tangent vectors to the parametric curves ( 8 = constant and cp = constant) are given by 5 - A A + + c cp =

-

b cos(8)sin(i)ex t b cos(û)cos(cp)~ Y ( 7 )

and the unit outward normal vector (for cos(fi)#o) by

n

- B - B

n = c/b

For cos(8)=0 the tangent vectors are dependent. In the software package this dependency is always checked and if necessary an error message is given. Points for which cos(8)=o should thus in principle be excluded. We may choose the movements, however, in such a way that these points never will become contact points. Therefore we will pay no further attention to this subject.

1.2 The rotation tensor.

The rotation tensor is defined as

+ T+

R = c e . ( 2 . 7 )

For -the -matrix repxesentation

-E

o f R- with respect- to the veetor base $- holas

I

from which follows

T

ox with the orthogonality relation & = I : -, T-,

E = & e.

..

5

From the orthogonality of 8 it fallows that only three of its nine

choosen as independent variables. The total rotation, leading from

e

5 to

:,

-.

may then be considered as the result o f 3 successive rotations, governed by the Cardan angles q ~ , 3, and w (see Fig. 4). The following relations apply:

5 5 5 5 5 c 5 5 yielding

Y

Fig. 4

with _P

where crp and srp stand for cos(rp3 and sin(rp) etc. With equations (16) and (?Sa3 we get for the rotation tensor

- b + e - +

R = (cw.cy, t sw.s$.sw)exext c$.srp.e e t X Y

- 8 + + +

t(sw.crp - srpp.s$.cw)e e t (-cw.sip

+

sw.s$.crp)e e

+

x z Y X

+ + + + - b +

tc$.cipe e t (-sw.sip

-

cw.s$.cip)e e - sw.c$e z x e t

Y Y Y Z

In- the present problem only motions kn the- xz-plane -occur In this case the rotation tensor is given by

i. e. ip -=- $ -= O.

- b + - b + + - + d + - b +

R = cos(w)(e e t e e

+

e e t sin(w)(exez-ezex) x x z z Y Y

with w the rotation about the y-axis as indicated in Fig. 4.

1.3 Formulation of the problem.

In the configuration described in the introduction (see also Fig.1) a constant external force fe in the negative z-direction is exerted on the center OA of the sphere. Furthermore, the translation of OA in the z- direction is prescribed. The following unknown quantities remain :

- the position of the contact points

-

the angle w

- t.he indeterminate multiplier D that is related to the load, necessary for

the prescription of the translation in z-direction.

With the assumption that - because of symmetry - only motions in the xz- plane will occur, the above quantities may be determined analytically with the use of the theory presented in WE 83-05, chapter 2.

This is the subject of section 1.4.

1.4 Analytical solution t o the Problem.

Because of the symmetry of the configuration at the contact points will apply cp = O on the cylinder and q~ = O on the sphere. The position of the contact points on these bodies is defined by h (see Fig.1) and 6.

1.4.1 Contact conditions.

____--_---

With ip = O it follows from equation [I) and from Fig.1 that

-P 4 -B

c =

sext

hez

and from equation (5) that

-B 4

n = - e .

X

119,

In a similar manner with cp = O it is derived from equations ( 6 ) and ( 9 ) that

. %

4 4

c = b c o ~ ( 6 ) ~ ~ t b sin(l)ez

and

If the prescribed translation of OA in positive z-direction is denoted by $,

then the translation vector

2

from OR to OA can be written as :

-+ 4 4

The contact condition (2.47) : "the contact points on sphere and cylinder should coincide" or

+ + +

c = a t R.c

yields with equations ( 1 8 ) , (191, (21) and (23) :

resulting in

b[ -1tcos (w) cos ( 8) +sin(") sin( 8

)]ex+

[ p-h+bIcos (")sin(

e)

-sin( w) cos ( 8)

1

]GZ=6.

A

o

(25)

The contact condition (2.48): "the unit normal vectors in the contact points have opposite directions" or

(261

t t

n =

-

R.n

results with equations ( I S ) , (20) and (22) in

-0

[ -Itcos ( w ) cos ( ;)+sin( w ) sin(

i)

]GX+

[cos ( w ) sin(

8)

-sin( w) cos ( 8)

]ez=d.

(27

1

J i I I

The constraint equations (25) and (27) are satisfied if the terms between

square bracket equal zero. Frcm and fOllOWS

A

*

sin(w) = sin(@) 4 w = 8 2 kii ( k = û , 1 ,

. . . .

1 .

1 . 4 . 2 The external loads on the sphere.

____-____-______--___---_---

The external loads acting on the sphere are given by

4

Ze

=

-

feez

( 2 9 )

The loads acting on the sphere due t o point contact are given by ( 2 . 8 7 ) and

( 2 . 8 8 ) according to A .+ 4 pn -B = - P x 4 f - -p R.n = k- 0 . 0 . d 9 - 9

and-because- it applies that c*n=o-it fel-luws- that +

$=

o . (34)

1.4.4 The loads on the sphere due t o kinematic constraints.

_________-__________---_-_---_-

The kinematic constraint resulting from prescribing the translation $ o f ûA

can be formulated as

- $ +

a.e = $

2

or according to (2.59) in the form k . ( a . R ) = O : 3

+ 4

k = a.e _z - $ = O .

From the requirement 6k=0 it follows that

and with (2.59) : k(&6a,R+öR) = O. Applying (2.60) to ( 3 6 ) yields + + a (at6a).ez- $ = O and, consequently

6Z.Z

_z =

o.

Relation (39) is the constraint equation o f the form

for the present problem. Obviously we have :

(2.60)

(39)

(2.61)

The loads due to the kinematic constraint become with the relations (2.112) and (2.113):

+

m = - I

0 . 3

-

6

r 3 I j -

1.4.5 Loads on the sphere due t o the line elements.

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ - _ _ _ _ _ - - - _ _ _ - - - _ - - - - _ - -

The constitutive relation for the line elements is chosen to be :

2 f(L) = A.L t EL

f : magnitude of the force in the element

L : actual length o€ the element A,B : constants.

The position vectors of the attachment points o f the line elements in the reference configuration are given by :

.+ + -+ s l = a x + slez ; + s = -bz 1 Y -b + - > + -4 s2= ae t s e (=sl) ; s2= bg

x

1 2 Y -> s3=

at

t s

e

x 3 % * ; -4 s3= -bgX

The force exerted on the sphere by element

where

2.

represents the unit vector directed along the line element and pointing from the cylinder to the sphere and

1

~ ~ ~ ~~~~~~~~~ ~

+ -b -4 + = ' ?

li = Livi = a t R . s . - s

1 i'

W i t h relations (18)' (23) and (44) to ( 4 6 ) inclusive €or the individual elements results :

element

1

:

(49)

-4 - # d

f l =

a t R . S ~ - sl=

-

bex- b e Y t ( @ - sl)e2 with R . g l = - be Y

2 2! 2b i- ( @

-

s l ) element 2 : - b + - b . + -4 + -9

1

= a t R.s

-

bex+ be t ( @

-

sl)gz with R.s2= be Y 2- "2= Y 2 ( 5 1 )

1 L2= L element 3 : -# + - # + f3= a f R.s3- s 3 =

-

b{l f cos(w))e _X t {b sin(w) t @ - s ₃

)e

₂ ( 5 3 ) A -4 -+ -+

with R.s3= - b cos(w)ext b sin(w)ez

2 2 2

Is = $b i 1 t cos(w)) t Ib sin(w)

+

$ - s3}

.'

(54)

The resulting force exerted by the line elements on the sphere is with the relations ( 4 7 ) and ( 4 9 ) to ( 5 4 ) inclusive found as :

3 3 + f ₁ = . E $.= -.E (A t BLi)fi 1=1 1 i=1 (55) b[2(A t % , ) + ( A 9 BL3)f1 i- COS(W))]~~- 4

- [ 2 ( A t BLl)(@

-

sl) -i- (A t BL3)jb sin(w) t

-

s3}]ez.

(56) -Ti% torque diie tö the line elements I s givSpi bY :-

3

-+

m 1 = i=l E

[R.S~)*Z~J

-+

= ( A t F3Ll)[- bey* {bext be

-

( 6

-

s1)~2t t bg

*

ib& - be

Y Y X Y

-e 3 -+

- !$ - s13e2 } ]

+

b(A t EL3)[{- cos(w)ex9 sin(w)ezl*!b(l t c 0 s ~ w 1 ) ~ ~

+ + * + 4

f

+

f l t fr+ fk = o +

e

-

feedz+ b[2(A t BL,) t (A t BL3){1 t c o s ( ~ ) f ] ~ ~ t

-

[2(A t BL,)($

-

s l ) t ( A t BL3)Ib sin(w)

+

@ - s3}lez- oe -9 -+

z

- pe - X e= o +

and the moment equilibrium gives

The component equation o f ( 5 7 ) and 158) are three equations f o r the unknown quantities p, w and u, yielding :

With h = 8 , 9 = w i 2kfi and a known value of @ the solution to the problem is known.

An araalysis Was done using the numerical values :

a = 2 cm, b =

1

cm, A =

1

N/cm, B = 2 N/cm

,

s l = 5 cm, s 3 = 6 cm and fe= 50 N.

2

In the Figures 4 . 4

-

4 . 6 of WE 83-05 the analytical and numerical results can be seen to be in excellent agreement.

PART I1

The following contains the listings of the user written subroutines that are needed for the numerical treatment of the problem, followed by a listing of the standard as well as the problem dependent input data (section 11.21.

Beforehand in section I I . l attention is paid to the calculation of áp, 6$

and 6w which are needed to update the rotation matrix in the subroutine UPDR.

11.1 Determination o f 6 ~ , 64 and 6w.

The axial vector 6 t of 6R .RC in general will. be a function of 6p, S$ and 6 w .

The components of

d

(with respect to

e)

are calculated in the program. Consequently 61p, 6$ and 6w may be formulated in terms o f the known components of

6.

A straightforward but rather tedious method to obtain the values of 6p, S$ and 6w is : determine the expression for Sg g in terms o€ p, $, w and their variations. Then use the knowledge that

,.,

T

which results in three equations that can be solved for ö q , ti$ and 6w

T

Another method circumvents the calculation of 6g g following :

From equation ( 1 6 ) we have

and comprises the

c p * w R = R .R .R and thus ( 6 7 ) 6 R = 'R .R*.RWö(p

+

Rcp.R* .RW6*

+

R 9 4 J w .R .R 6w r c p

,*

,w

where the subscript i a denotes differentiation with respect to a. From (669

and ( 6 7 ) we find :

C

with the scew-symmetric tensors A"= Raa. Ra having axial vectors

afa.

Using the right hand side of equation (68) we will elaborate the expression

I

c +

6R.R .v = 6if.f for arbitrary

c.

It follows that ( 7 0 ) cp + +P*+ A .v = 6li v C ( 7 1 1 c p * C R@.A*.Rv .v" = R .A .(RQ ,v") = R Y { ~ ~ * * ( R Q ~ . ; ) ~

Because a rotation tensor R and its adjugate are identical, i.e. Ra.($

*

₌ R.(;

*

4)

= (R.;)

*

(R.;)

we get

In a similar way we Eind

Combining (68), ( 6 9 1 , (701, (73) and (74) yields eventually (in matrix representation):

where the components o f the axial vectors 6fla follow from

11.2. LISTINGS OF USER WRITTEN SUBROUTINES SUBROUTINE EXTBEL(N1NC)

$INSERT WFW>FONS>KNIE.SOURCE.DIR>HEADER C

C FORCE: -FE=-USER(250) IN Z-DIRECTION C FE(lf3)=-USER(250) RETURN END C C C $INSERT W ~ W > ~ O N ~ > K N I E . S O U R C E . D I R > H E A D E R C C C SUBROUTINE FITSUR

PAROP(l)=A AND PAROPD(l)=B ARE STORED IN RESTARTFILE PAROP(l)=USER(200) PAXOPD(l)=USER(201) RETURN END C C C SUBROUTINE GEMET(NOPP,PKP) $INSERT W F W > ~ O N S > K N I E . ~ O ~ R C E . D I ~ > H E A D E R IF (NOPP.EQ.2) GOT0 1000 C

c

C C C C C

c

C C C C C C C C C C C

DETERMINATION OF GEON~TRIC Q~7ANTIT~ES FOR BODY 1 (CYLINDER) AT CONTACT

POINT IKP: ~ ~~~~ ~

1). COMPONENTS OF POSITION VECTOR OF CONTACT POINT ARE STORED IN C(IKP,1-3)

2). COMPONENTS OF TAMGENT VECTORS IN CONTACT POINT ARE STORED IN

3). COMPONENTS OF UNIT NORMAL VECTOR ARE STORED IN N(IKP,1-3)

4 ) . CO~PONENTS OF RECIPROCAL BASE VECTORS ARE DE~ERMINED BY MEANS OF THE

5). COMPONENTS OF CURVATURE MATRIX ARE STORED IN K(IKP,1-3,1-3). CYLINDER = BODY 1; U=U(l,I] =,9,,V=U(lr2)=2

SPHERE = BODY 2: U=UD(l,l) = 4 , V=UD(1,2)=8 CU(IRP,1-3) AND CVfIKPrI-3:.

SUBROUTINE INVBAS AND STORED IN BU(IKPr1-3) AND BV(IKP,1-3)

FOR THE SURFACE PARANETERS U, V, U AND V APPLIES:

BODY 1 = CYLINDER

1). CALCULATION OF THE COMPONENTS OF THE PûSITION VECTOR OF CONTACT POINT 1 (THE ONLY CONTACT POINT) ACCORDING TO EQUATION (1):

COSF=DCOS(U(l,l)) SINF=DSINIU(l,l)) RA=PAROP(I)

C(l,l)=RA*COSF C(1,2)=RA*SINF

C C C C C C C C C C C C C C

c

C C C C C C C C C C

2). CALCULATION OF THE COMPONENTS OF THE TANGENT VECTORS ACCORDING TO EQUATION (2) :

CU(I,A)=-RA*SINF CU(I,2)=RAXCOSF CU(lf3)=0.D0

AND ACCORDING TO EQUATION ( 3 ) : CV(I,l)=O.DO

CV(If2)=O.DO CV(i13)=1.D0

3). CALCULATION OF THE COMPONENTS OF THE UNIT NORMWE VECTOR ACCORDSMG TO EQUATION ( 5 ) :

N(I,l)=-CûSF M(lf2)=-SINF N(If3)=O.DO

4 ) . CALCULATION OF THE COMPONENTS OF THE RECIPROCAL BASE VECTORS:

CALL INVBAS(NOPP,IKP)

IF (.NOT. LCHECK(2,IKP)) GOT0 10000

REMARK: LCHECK (2,IKP)zTRUE IF THE LOCAL VECTOR RASE IN IKP IS REGULAR.

5 ) . CALCULATION OF THE COMPON~NTS ~ OF ~~ THE ~ ~ ~ R V ~ T U ~ E ~ ~~~ MJATR ~

~~~ ~~~~ ~~

~~~ ~~ ~

THE C ~ T g N p 3 ~IS+GJYEtJ ~Ex: A ~ ~ ~ ~

K = R

B

B

+E

B

B +K

B

B +K

BvP,

AND 4

IN'

-M~~ P I

X'#EYRE

s

EEYAY

I ~ N

5 . "

-y3+--qv

+f" VV,

K = K e e e p @ e * e - K uv, p p T z ~ us v uv-uv R

+ 3cv -+ 9

K = e o K o e = R +R +IJ +IJ WITH

KuV

= K vee@Lipvee -uv uv, z ..usv5 ,

E I S THE MATRIX REPRESENTATION OF THE DYAD WITH

-UV

%v IS CALCULATED BY THE SUBROU5$IHE DYADPR.

- n * c = n T c ~URTHERMORE WE HAVE Kuv = n o

-

-

+T T+ AHD 6 = e 63 = f3 e FOLLOWS: - s .l... +

-#

T 4 +T"- uv 5 ,uv' au CONTACT IX

E.

~ KNQWH B, 4 POINT

WITH KNOWN n AND c THIS SCALAR PRODUCT MAY BE CALCULATED WITH THE SUBROUTINE RINWK. s *uv :OR TH$ CYLINDER APPLfES: -b = -a cos(cp)e

-

a sin(cp)e 4 K U u # O.

3

'

TsC9~oRED IN KLADVT: Y z uu C

C C C C C C C + +

aCZ

(= = c = - -

-

6

+ Kuv = O. uv yrz _$$ c = c

=--%.a

+ K V U =

o.

vu Zlp

9;

c = c = - - ' - - d + K v = O . F ~ E T H E ~ ~ R E S ~ ~ ~ T PROBLEM

EESULTS

: + -b .+ -b T T U U uu

..

,UU "U" u K = K

CALCULATE AND STORE IN KLADNI:

WITH KUU = n c _AND

_Buu

_{= $ fl}

_.

-

C CALL DYADPR(BU~BU~KLADMI~IKP,IKP,I,5,5,1) C C C C C C C C C C 1 000 C C C C C C C C C

CALCULATE KuU AND STORE IN KLADV2(1,1):

CALL KINWK(N,KLADV1,KLABV2,IKP,I,1,1,1,1,1) CALCULATE = KUUBUU: SCA=KLADV2(1,1) CALL SCA~(KLADNI,SCA,K,liIKP,1,5) GOT0 10000 BODY 2 = SPHERE

1). POSITION VECTOR ACCORDING TO EQUATION (6):

COSFD=DCOS(UD(I,I)) SINFD=DSIN(UD(l,l)) ~~ COSTD=DCOS (üD ( î

' 2 )

1

RB=PAROPD(I) CD(I,I)=RB*COSTD*COSFD CD(lr2)=RB*COSTD*SINFD CD(I,3)=RB*CINTD ~~~~~~~~ ~ ~~ STNTD=ESTN(UE!l,2!)

2). TANGENT VECTOR ACCORDING TO EQUATION (7):

CUD(I,l)=-RB*COSTD*SINFD €UD(1,2)=RB*COSTD*COSFD CUD(1 I 3)=O.D0

AND ACCORDING TO EQUATION ( 8 ) : CVD(1,1)=-RE*SINTD*COCFD

CVD(l12)=-RB*CINTD*SINFB CVD(II3)=RB*COSTD

3). UNIT NORMAL VECTOR ACCORDING TO EQUATION (9):

1100

DO 1100 I=1,3 ND!Z,I!=€D[I,I)/RB CONTINUE

C C

C

4). RECIPROCAL BASE VECTORS:

C C C C C C C C C C C C C C C C C C C C C CALL INVBAS(NOPP,IKP)

IF (.NOT. LCHECK(2,IKP)) GOT0 10000

5 ) . C U ~ V A T ~ R ~ ~ M A T R I X : 3 1 c = c 3 c = c = - - 3 1 -b -t c = c THIS ,YIELDS: -b ac

= -A = -b cos8 cos@ ex - b cos8 s i n @

z

1 + 3 = c Y

*

- b s i n 8 s i n @ e

-

b s i n 8 cos@

e

uu P P

%i!

uv tp8 X Y vu = - - - -b cos8 C O S ~ e x - b cos8 s i n @ e - b s i n 8

eZ

vv 88

ae

Y CALL MNULS(KLADV1,1,1,1) KLADV1(I,l)=-RB*COSTD*COSFD KLADV1(I,%)=-RB*COSTD*SINFD KLADVl(I13)=0.D0 *

STORE Kuu IN KLADV2(1,1):

CALL KINWK(ND,KLADVl,KLADVS,IKP,1,1,1,S,1,1)

A

CALCULATE

-uu

AND STORE IN KLADMl: IN KLAûM4{lf..,..): STORE %U = KUU %u

SCA-KLADVS(1,l)

CALL SCAMPI (MLADMI

,

SCA, KLADMI4

1

,

I

1

12 1 CTORE

cuv

IN KLADV1:

5

CALL NNULSiKEAE)ViI1,l,l) KLADVl(l,l)=RB*SINTD*SINFD KLADV1(1,2)=-RB*SINTD*COSFD KLADV1(Ir3)=O.DO

STORE KUv IN KLADV2(1,1):

CALL KINWK(ND,KLADV1,KLADV.2rIKP~lfl,lf51111~

STORE

kv

IN KLADMl:

CALL D Y A D P R ( B U D f B V D , K L A D M l ~ I K P ~ I K P ~ l f ~ ~ 5 ~ l )

STORE

%v

= KUv%, IN KLADM4(2,..,..):

C C C C C C C C C C C C C C C C C C C C C C C ~~ ~ 10000 SCA=KLADV2(1,1) CALL SCAMM(KLADM1,SCA,KLADM4,1,2,1,12) STORE

yu

IN KLADMI: CALL DYADPR(BVD,BUD,KLADM111KP,IKPII15fS,I) A A A A A

STORE

"vu

= K B = K B IN KLADM4(3,..,..): vu-vu U V T U CALL SCAMM(KLADM1,SCAIKLADM4,1,3,1,12) STORE

tvv

_..

IN KLADVI: CALL MNUL2(MLADVI,1,1,1) KLADVI(1,?)=-RB*COSTD*COSFD KLADV1(1,2)=-RB*COSTD*SINFD KLADV1(II3)=-RB*SINTD STORE Kvv IN KLADV2(1,1): CALL KINWK(ND,KLADV1,KLADV2,1KP~lIlI~z5~l,l) STORE B IN KLADMI:

CALL Di%PR(BVD

,

BVD I KLADMI I IKP IKP

,

1 S,5

,

1

1

STORE K = K B IN KLADM4(4,..,..): S C A = K L a b ( 1

,hiYTv

CALL S C A ~ ( K ~ A D M 1 , S C A , ~ L A D M 4 , 1 , 4 , 1 , l 2 ) * * A CALL ~ ~ L U S M [ K L A D M 4 , K L A D M 4 ~ K L A D M 4 1 1 1 2 , 5 , 1 2 f 1 2 f I ~ ~ CALL MPLUSM(KLADM4,KLADM4,KLADM4,S13,6,12~12,12~ CALL M~LUSM~KLADM4,KLADM4~KD,6,4,1,12,12,S~ GOT0 10000 RETURN END C

c

C C SUBROUTINE INUSER $INSERT WFW>FONS>KNIE.SOURCE.DIR>HEADER C C

1

O0 C C FORMAT(A2)

C RESPECTIVELY C READ(IFINF,IOO)BEGIM READ(IFINF,*) USER(2001 READ(IFINF,*) USER(201) READ(IFINF,lOO)EIND C

C STORE VALUES OF PRESCRIBED TRANSLATIONS í3 IN USER(225:224+MAXINC)

c

READ(IFINF,100)BEGIN MAXINC=IUSER(13) DO 90 I=I,MAXINC READ(IFINF,*) USER(224iI) READ(I~IN~,IO~)EIND 90 CONTINUE C C STORE FE IN USER(250) C REA~(IFINF~IO0)~EGIN READ(IFINF,*) USER(250) READ(IFINF,IOO)EIND C

C STORE STIFFNESS CONSTANTS A AND B OF THE LINE ELEMENTS IN ARRAY USER C ACCORDING TO:

C A(i=1) R(i=2)

C line elem. nr. 1 user(260) user(261)

I' 3 user(264) User (265 ) C C It 2 user(262) user (263 1 I 1 I< 11 II READ~I~INF~IOO)BEGII~ DO 110 ILYN=l~~AXLYN READ(PFINF,*) (USER(259+2*(1LYN-1)~1~,1=1~2) R E A D ( I F I N ~ , I O O ~ E ~ N ~ RETURM END ~ ~ ~~ ~ ~ ~ ~~~ ~ ~~ ~~ ~~ ~ ~ M ~ ~ ~ Y N = I U S E R ~ ~ ~ ~ ~ ~ ~ 110 CONTINUE C C C

e

$INSERT ~ F ~ > F O N ~ > ~ N I E . S O U R C E . ~ I ~ > H E A D E R C SUBROUTINE KINECO(IKICO,NSNC)

WE HAVE GOT ONLY ONE KINEMATIC CONSTRAINT:

a t - $ +

a e e = $ a R $ z a t e - $ = O

cx, * 6a

+

Ti, Biï = O A -L

z +Z +

6kl = O + (a+6a]te -@-(see -6) = O

6k =

5)

M$Y ALS9 PE WRITTkN AS

+ I

I$ THE FJRES~NT

CAS^

WE HAVE

úl

= éz

6 q , = 6 + D 3 = D 4 = 0

OTT1 =

6

+ D, = D, = O

IF (IKICO.NE.1) GOT0 1000

= o

LCHECK(S,I)=.TRUE. DO 100 IKICO=2,5 LCHECK(S,IKICO)=.FALSE. ALPHA( 1,3)=1 .DO END 100 CONTINUE 1000 RETURN C C SUBROUTINE REUSER $INSERT WFW>FONS>KNIE.SOURCE.DIR>HEADER C

C ROUTINE TO READ SOME PARAMETERS FROM THE C PARAMETERS SHOULD BE READ IN ARRAYS USER C AND 25 ON, RESPECTIVELY

C PARAMETERS SBOUED BE READ FROM FILE WITH C HAS BEEN INITIALIZED IN ROUTINE INIT2

C

INPUTFILE; IN CASE OF RESTART AND IUSER FROM POSITIONS 200

FORTRAN U ~ ~ T N U M B E ~ IFINF WHICH

C

C WRITTEN TO THE RESTARTFILE C C

c

C C C

PARAMETERS TO BE READ: USUALLY THOSE WHICH ARE NOT AUTOMATICALLY

IN THIS SIMPLE MODEL WE WILL NOT USE THIS ROUTINE ( WE WON’T RESTART )

RETURN END C

c

SUBROÜTÏNZ SPLÏTC-( NOPP I IKP

1

$INSERT W~W>FONS>KNIE.SOURCE~DIR>HEADER

C

C

~ ~~~~

~

C ~ETERMINATIO~ OF THE SURFACE C~ORDINATES

u

AND

v

FROM A GIVEN POSITION

c

VECTOR. THIS POSITION VECTOR BELONGS TO THE SPACE-FIXED BODY IF NOPP = 1 C

C TACT POINT WE CONSIDER. THE POSITION VECTOR I S GIVEN IN ARRAY C IN POSI- AND TO THE MOVING BODY IF NOPP I=. 2. THE POINTER IKP TELLS US WHICH CON-

C C C C C C

c

C C C C C C C C

TIONS C(IKP,1..3) IF NOPP = 1, AND IN ARRAY CD IM POSITIONS CD(IKP,1..3) IF NOPP = 2. THE SURFACE COORDINATES MUST BE PLACED IN ARRAY U IN POSI- TIONS U(IKP,1..2) IF NOPP=I AND IN ARRAY UB IN POSITIONS UB(IKP,1..2) IF NOPP = 2

THIS ROUTINE IS ONLY NEEDED IF:

- THE LENGTH OF A LINE ELEMENT BECOMES ZERO. THEN IN ROUTINE CHELYN A CONTACT POINT IS GENERATED. THE INSERTION POINTS OF THE LINE- ELEMENT

THEN BECOME THE COORDINATES OF THE NEW CONTACTPOINT, RUT THEY MUST BE SPECIFIED IN TERMS OF THE SURFACE COORDINATES

- _{IN THE INPUT DATA A CONTACT POINT IS GIVEN IN ITS COORDINATES AND C}

C C C EMPTY BODY C C C

C IN OIJR PROBLEM THIS ROUTINE ISN'T NEEDED C

C

C

C

IF THESE CIRCUMSTANCES DON'T OCCUR THEN THIS ROUTINE CAN HAVE AN

C C RETURN END

c

C SUBROUTINE SPLITR $INSERT WFW>FON~>KNIE.SOURCE.DIR>HEADER C C

C ROUTINE TO DETERMINE THE ROTATION PARAMETERS FROM THE ROTATIONMATRIX R

C

C SECOND INDEX IS A POINTER TO THE ROWS OF THE MATRIX

C THE ROTATION PARAMETERS MUST BE STORED IN ARRAY PHI IN POSITIONS

THIS MATRIX IS GIVEN IN ARRAY R IN POSITIONS R(1,l

. . .

3 , 1

. . .

3 ) WHERE THE

C PHI(1,1..3J

C THIS ROUTINE MUST BE GIVEN IF WE DIDN'T

C IN THE INPUTFILE BUT THE RO~ATIONMATRIX C HAVE AN EMPTY BODY

C

c

C

C IN OUR MODEL WE WON'T NEED THIS ROUTINE

C C c ~ ~ ~~ RETUR~ END C C SUBROUTINE UITINC~NINC'ITE~9 $INSERT WFW>FONS~KNIE.SOURCE.~IR>HEADER $INSERT SYSCOM>A$KEYS C C C C

GIVE THE ROTATIONPARA~ETERS ITSELF,ELSE THIS ROUTINE CAN

C ROUTINE TO GIVE SOME OUTPUT IN A DISKFILE FOR THIS INCREMENT C OUTPUT MUST BE WRITTEN TO FORTRAN LOGICAL FILENUMBER IFUITF C

C VARIABLES

C

C

C NINC = CURRENT INCREMENTNUMBER

ITER = NUMBER OF ITERATIONS IN THIS INCREMENT 10 20 1 O0 110 120 130 140 150 C C C START OUTPUT

__-___----__---

WRITE SOME RELEVANT PARAMETERS

WRIT€(IFUITF,I)

THIS FORMAT PROVIDES FOR A NEW PAGE AT THE PRIME PRINTER FOR EACH INCREMENT

FORMAT(2H ,/,IHlI WRITE(IFUITF,IO) FORMAT ( /

,

1 X i 35 ( 1 H$ ) )

W R I T ~ ~ ~ ~ ~ ~ T ~ ~ 2 0 ) ~ ~ ~ ~ FORMAT(/,IX,26HUITVOER VAN I~CR€MENT NR :,12)

WRITE ( IFUITF,

I

o)

WRITE(IFUITF,lOO) FORMAT(//rlX,19H~LGEMENE GROOTHEDEN) WRITE(IFUITF,l10) FORMAT(lX,lS(lH=)) WRITE(~FWITF~l2O)NINC

FORMAT(lX,l~HINCR~MENT NUMMER :,I21 WRITE(IFUITF,130)ITER

FORMA~(1X~I~HAANTAL ITERATIES :,12) WRITE(IFUITFf140)NKONP

FOR~AT(l~rl~HAANTAL KONTAKTEN :f12) WRITE(IFUITFf~5~)NKINC0

FORMAT(1X~l~HAANTAL KOMDITIES :,I21

WRITE PROCESSORTIME FOR THIS I N C R ~ ~ ~ N T

~~ ~ C C 160 170 C C C C 200 210 220 TTYD=DBLE(FLOAT(ITOTYDf)/lOO.DO T Y D I ~ C = ~ B ~ E ( F L O ~ T ~ I ~ C T Y D ) ) / ~ O O . D O WRITEIIFUITFfi60)TYDIN@

FORMAT(lX, 28HCPU TYD VOOR DIT INCREMENT :,F11.4,2Xr4HSEC.) ~ R I ~ E ( I F U I T F r ~ ~ O ) T T Y D

FORMAT(lX, IGHTOTALE CPU TYD :,F11.4,2X84HSEC.)

WRITE GEOMETRICAL PARAMETERS

WRITE(IFUITF,200) FORMAT(//,IX,23HGEOMETRISCHE GROOTHEDEN) WRITE(IFUITFf210) FORMAT(1Xf23(1H=)) w R r r E ( m n T F , z ~ ó j FORMAT(/,lX,SHLICHAAM 1)

230 260 245 250 270 280 2 90 5 O0 600 700 C ë ~ C (3 C i O00 I010 1020 1030 1031 1040 1 O50 1060 1070 WRITE(IFUITF,230) FORMAT(IXf9(1H-)) DO 500 IKP=1,5

IF ( .NOT. LCHECK( 1

,

IKP) ) GOTO 500 WRITE(IFUITFf26O)PKP

FORMAT(/,IX,I6HCONTACTPUNT NR :,12) W R I T E ( I F U I T F , 2 4 5 ) ~ U ( I K P , I ) , I = 1 , 2 )

FORMAT(/,1Xr24HOPP. PARAMETERS U EN V :,T30r3(F11.4,1X)) WRITE(IFUITF,250)(C(IKP~I~,I=lr3~

FORMAT(1Xf29HPOSITIEVECTOR KONTAKTPUNT C :,T3Of3(F11.4,1X) 1 WRITE(IFUITF,270)

FORMAT(IX,26X,lH")

WRITE(IFUITF,280f(N(IRP,I)~I=1,3)

FORMAï(1~,28H~ORMAAL OP HET OPPERVLAK N :,T30,3(F11.4,1X)) WRITE(PFUITF,290) FORMAT(1Xr25X,1H") CONTINUE WRITE(IFUITF,GOO) FORMAT(//,IX,9~LICHAAM 2) WRITE(IFUITF,230) DO 700 IKP=I,5

IF ( .ROT. LCHECK(1,IRP) ) GOTO 700

WRITE(IFUITF,26O)~K~ W R I T E ~ ~ ~ U ~ T F r 2 4 5 ) ( U D ( I K P f I ) í I = l f 2 ~ WRIT~(IFUITF,25O~(CD(IKP,I),I=l,3) WRITE(IFUITF,270) W R I T E ( I F U I T F , 2 8 0 ~ ( N D ( I K P , I ~ , I = 1 , 3 ~ WRITE(IFUITF~290) C O ~ ~ I N U E ~~~~~~ ~ ~~~~~~~~~~

WRITE KINEMATICAL PARAMET~RS

WRITE(I~UITF~IOO0) F G ~ ~ ~ T ~ ~ ~ ~ l X ~ 2 3 H K I N E M A ~ ~ S ~ ~ E ~ R O O ï H E D ~ ~ ) WRITE(IFUITF,210) W R I T E ( I F U I T F , 1 0 1 0 ) ( A ~ l f I ~ , I = l ~ 3 ~ FOR~T(/llX,20HTRANSLAT1EVECTOR A :,T30,3(F11.4,1X)) WRITE(IFUITF,1020)(R(l,lfJ~rJ=l,3~ FO~MAT(lX,l7HROT~TIEMATRIX R :,T30,3(F11.4,1X)) WR1TE(1FUITF,1030)[R(1,2,J),J=1,3) FORMAT(1Xr14X,1H-,T30r3(Fll.4,1X)) WRITE(IFUITF,1031)(R(lí3,J~fJ=1,3) FORMAT(16X,T30,3{Fll.4,'lX)1 WWITE(IFUITF,1040)(PHI(lfI~,I=l13~ FOR~T(lX123HROTATIEPARAMETERS PHI :,T30,3(F11.4,1X)) WRITE(IFUITF,l050) FORMAT(lX, 19X, IH") WRITE(IFUITF,1060)ROTSAS

FORMAT(IX,S4HROTATIE ROND SCHROEFAS :,T30fF11.4) WRITE(IFUITF,1070)(SAS(I),I=lf3~

FORMAT(IX,20HRICHTING SCHROEFAS :,T30,3(Fl1.4,lX))

1080 1 O90 C C C C C 1200 1210 1220 1230 1240 1250 1260 1300 1310 1320 ~~ ~~ 1350 1360 C

c

C C C 1400 1410 1420 1430 1440 1450 1500

FORMAT(lX~28HTRANSLATIE LANGS SCHROEFAS :,T30fF11.4) WRITE(IFUITF,lO9O)(STVSAS(I),I=l,3~

FORMAT(1Xr23HSTEUNVECTOR SCHROEFAS :,T30,3(F11.4,1X))

WRITE LOADS DUE TO CONTACT POINTS

WRITE(IFUITF,1200) FORMAT(//,lXf20HBEL~STINGSGROOTHEDEN) WRITE(IFUITF,1210) FORMATflX,20(1H=)) WRITE(IFUITF,1220) FOR~T(/,lXi15HKONTAKTKRACHTEN) WRITE(IFUITF,1230) FûR1VIWT(IX,I4(iH-)]

IF ( NKONP .EQ. O ) GOTO 1350 DO 1260 IKP=1,5

IF ( .NOT. LCHECK(1,IKP)) GOTO WRIT~(IFUITF,l24O)IKP,PK(IK~) 2 ~ O R ~ T ( l X , ~ l H C O N T A ~ T ~ U N T f 1 2 ~ 2 X , 4 H p :,T30fF11.4) WRITE(IFUITF,1250) FORMAT(1Xf16X,lHk) CONTINUE WR1TE(1FUITF,1300)(FK(1,1~,1=1,3) FOR~T(/,lX,25HTOTALE KONTAKTKRACHT F :,T30,3(Fl1.4,lX)) WRITE(IFUITF,1310) FORMAT(1X,21X,2HNk) WRITE(IFUITF,1320)(MK(l~I)~I=l,3) ~ O R M A T ~ ~ ~ , 2 5 ~ T ~ T A ~ ~ CONTACT MOM€^^ M :,T30,3(F11.4,1X)) WRÏTTIIFUITF, 1 3 1 3 ) ~ GOTO 1400 WRITE(IFU~TF,13~0) FORNAT(/,lX,2OHGEEN CONTACTKRACHTEN) ~~ ~~ ~~ ~ ~~~~ ~ ~ ~ ~ ~ ~~~~~~

WRITE LOADS DUE TO KINEMATIC CONSTRAINTS

WRITE(IFU~TF,1410)

~ O R ~ A T ( / , ~ X f 3 ~ H B E L ~ S T I N ~ ~ N TGV KINEMATISCHE CONDITIES) WRITE(IFUITF,1420)

FORMAT(1Xr38(lH-))

IF ( NKINCO .EQ. O ) GOTO 1600 DO 1450 IKICO=1,5

fF ( .MOT. LCHECK(5,IKICO)) GOTO 1450 WRITE~IFUITF,l43O)IKICO,SIGMA(IKICO~ F O R M A ~ ~ l X ~ 2 l ~ K I N E N A T I S C H E CONDITIE,I2,2X,8HSIGMA :,T35,F11.4) WRITE(IFUITF,1440) FORMAT(lX,30X,lHj) CONTINUE WRITE(IFUbTF,?500)(FR(l~I),I=l,3) FOXMAT(/,fXf27HTOTALE RESTRAINTKRACHT F :,T30,3(F11.4,1X)) WRITE(IFUITF,1510)

WRITE LOADS FOR LINE-ELEMENTS 1510 FORMAT(1X,23Xl2H"r)

1520 FORMAT(IX,27HTOTALE RESTRAINTMOMENT M :,T3Of3(F11.4,IX)) WRITE(IFUITFf1520)(MR(l,I~,I=l,3~ WRITE(IFUITF,1510) GOT0 1700 1600 WRITE(IFUITF,1610) 1610 FORMAT(/,lXI25HGEEN RESTRAINTBELACTINGEN) C C C C

c

1700 WRITE(IFUITF,1710) 1710 FORMAT(//,1X129HBELASTINGEN TGV L I ~ N E L ~ M ~ N T E N ) WRITE(IFUITF,1720) 1720 FORMAT(lX,29(1H-)) MAXLYN=IUSER(l4) DO 1740 I=?,MAXLYN WRITE(IFUITF,1730)I,LVEER(I),FLYN(I) 1730 F O R ~ T ( l ~ , 7 H E L E M E N T I 1 3 , 2 X , l ~ H L E N G T E EN KRACHT,T30,2(FI1.4,1X)) 1740 CONTINUE

1750 FORMAT(/,I~,3OHTOTALE KRACHT TGV ELEMENTEN F ,T35,3(F11.4,IX))

1760 FORNAT(1X,28X,2HN1)

1770 FORMAT(IX,30HTOTALE MOMENT TGV ELEMENTEN M ,T35,3(F11.4,IXf) C

c

C W R m Y EXTERNAL -LoÄBs C C WRITE(IFUïTF,I900) 1900 ~ O R M A T ( / / , l X , 2 2 H ~ I T ~ E N D I G E BELASTINGEN) WRITE(IFUITF,lSIOI 1910 FORMAT(IX,20(1H-)) WRITE(IFUITF,1920)(FE(lfI),I=l~3) 1920 ~ 0 R M A T ( / , 1 X , 2 1 H U 1 T ~ E N D 1 G ~ KRACHT F : ,T30,3(F11.4,1X)) WRITE(IFUITFf1930) 1930 FORMAT(IX,18X,2H"e) WRITE(IFUPTF,1940)(ME(l~I),I=l,3} WRITE(IFUITF,1930) WRITE(IFUITF,1750)(FL(l,I),I=I,3) WRITE(IFUITF,1760) WRITE(IFUITF,1770)(ML(1,I),I=1,3) WRITE(IFUITF,1760) - - ~~ ~ 1940 FOR~T(lX,2lHUITWENDIGE MOMENT M :,T30,3(F11.4,1X)) C C C

c

C C C C RETURN END

C SUBROUTINE UPDA(N1NC) $INSERT WFW>FONS>KNIE.âOURCE.DIR>HEADER C C C A IS PRESCRIBED C

z

C

C A AND A ARE UPDATED BY MEANS OF ITERATIVE CORRECTION

C

x

Y C

c

C C A(lf3)=USER(224tNINC) DO IO0 I=l,2 A(1,I)=A(I,I)tDA(1,1) 100 CONTINUE RETURN END C C C C

SUBROUTINE VEREN ILYN,FVEER,STIFF $INSERT WFW>~ONS>KNIE.SOURCE.DIR~HEADER C C C ~

-c

C C C C C C C C C C C C C C C C C C C C C

c

ROUTINE TO CALCUHITE FORCE FVEER ANIT ELACTTC CTIFF%ESS STIFF OF LINE-ELE~~NT NUMBER ILYN FROM THE LENGTH OF THE LINE-ELEMENT GIVEN IN ARRAY LVEER IN POSITION LVEER(ILYN)

IN OUR EXAMPLE WE HAVE

FVEER A(L - L ) t B( L

-

L )**2

O

o

STIFF = A i- 2*B*(L - L )

O

FOR LINE-ELEMENT ILYN THE PARAMETERS B , B AND L ARE STORED IN ARRAY

IN POSITIONS USER(260+2*(ILYN-i)

....

262+2*(ILYN-I)) RESPECTIVELY

O

OMIT THIS IF LENGTH .LT. O

C C C C

I

O0 C 200 RLV=LVEER(ILYN) FVEER=RA*(RLV-RLO) t RB*(RLV-RLO)**2 STIFFzRA t 2.DOXRB*(RLV-RL0) GOT0 200 FVEER=O.DO STIFF=O.DO RETURN END

c

C C C C $INSERT W F W > F O N ~ > K N I E . S O U R C € . D I R > ~ ~ A D E R C

C UPDATING OF CARDAN ANGLES AND ROTATION MATRIX

C SUBROUTINE UPDR(NINC1 COSF=DCOS(PHI(l,l)) SINF=DSIN(PHI(I,I)) COSP=DCOS(PHI(1,2)) SINP=DSIN(PHI(l,2)) COSO=DCOS(PHI(1,3f) SINO=DSIN(PHI(1,3)) ~~ ~~~~p ~~~ ~~ ~~~ ~~ C

c

DETERMINATIOH OE THE CBAEsGES OF THE CARDAN ANGLES: pt%p: 6d, 6W :

C ~~ ~ ~ ~ p ~ ~ ~ ~ p ~ DPHI=SPNP*(COSF*DPI(lf2)tSINF*DPI(l,l))/COSP-DPI~~I~~ DPSI=COSF*DPI(l,l)-SINF*DPI(1,2) DOM=~COSF*DPI(1,2)+SINF*DPI(l~l))/COSP C C C

UPDATING OF CARDAN ANGLES

SCALING O F ANGLES BETWEEN O AND 2r. COSF=DCOS~PHI(Z,1)) SINF=DSINIPNI(lI1)) CALL CALCUH(COSF,SINF,PHIH) PHI(l,I)=PHIH COSP=DCOS(PHI(1,2)) SINP=DSIN(PHI(1,2)) CALL CALCUH(COSP,SINP,PSfH) PHI(1,2)=PSIH COCO=DCOS(PHI(1,3)) SINO=DSIN(PHI(1,3))

10 20 CALL CALCUH(COSO,SINO,OMEG) PHI(1,3)=OMEG IF(DABS(SINP).GE.USER(55)) GOTO 10 LCHECK(G,l)=.TRUE. CALL VULRMP GOTO 20 LCHECK(G,l)=.FALSE. RETURN END C C C C $INSERT WFW>FONS>KNIE.SO~RCE C C C C C SUBROUTINE VIILRMP DIR>HEADER

ROUTINE TO CALCULATE THE R O ~ A T I O N ~ T R I X FROM ITS ROTATIO~PARAMETERS PHI

w

THE MATRIX R IC STORED IN ARRAY R IN POSITIONS R(lfl..3$1..3) WITH THE SECOND INDEX AS A POINTER TO THE ROWS OF THE MATRIX

THE ROTATION PARAMETERS PHI ARE STORED IN ARRAY PHI 18 POSITIONS PHI(1,l. . 3 ) C C C C C C AMGEES C C

c

C

e

THE THREE ROTATION PARAMETERS ARE THE ANGLES PH1,PCI AND OMEGA

~

c-

SEE SECTION I,2 -TO GET TKE-PRECIFX F ~ R ~ U ~ ~ I O ~ OF IN TERMS OF THESE

C C COSP=DCOS(PHI(l,l)) SINF=DSIN(PHI(l,lf) COSP=DCOS(PHI11,2)) SINP=DSIN(PXIf~,2)) COSO=DCOS(PHI(1~3)) SINO=DSIN{PHI(I,3)) R ~ ~ , ~ , ~ ~ = C ~ S F " ~ O S O + S I N F " S I N P " S ~ N O R(1,1,2)=SINF*COSP R(1,1,3)=COSFkSINO-SINF*SINP*COS0 R(1,2,1)=COSF*SINP*SINO-SINF~COSO R(l,2,2)=COSF*COSP R(1,2,3)=-1.DO*(SINF*SINOtCOSF*SINP*COSO) R(1,3,1)=-l.DO*COSP*SINO K i 3

,

2 j =C IW R(1,3,3)=COSP"C0S0

C

RETURN

11.3. INPUT DATA

STANDARD INPUT DATA BEGIN O END REGINTRANS I.DO,O.DO,-I.DO ENDTRANS BEGINROT

1

O.DO,O.DO,O.DO ENDROT

BEGIN CONTACT POINTS

1

I

O.DO,-1.DO 1

O , DO, 0. DO

END CONTACT POINTS EEGININC 6 ENDINC BEGINSINGCP O. OIDO ENDSINGCP BEG1 NCONV O.OOIDO ENDCONV BEGIN1 TER 50 END ITER REGINSINGR O. 99DQ ENDSINGR BEGIN LINE CP - ~ ~~~- ~ ~ - ~ _ ~ _ ~~~ ~-~~~~ ~ ~~ ~- ~~~

0.01DO END LINE CP BEGIN LINE 3 2.DOrO.D0,5.DO 2.DO,O.D0,5.DO 2.DO,O.DO,6.D0 O.DO,-l.DO,O.DO O.DO,I.DO,Q.DO -l.DO,O.DO;O.DO END LINE

END OF STANDARD INPUT DATA

PROBLEM DEPENDENT INPUT DATA BEGIN RADIUS SPHERE AND CYLINDER 2 .DO

1

.no

END RADIUS SPHERE AND CYLINDER BEGINKINCO: PRESCRIBED TRANSL.

-1

.DO ~~ O.DO ~

1

.DO 2.DO 3 .DO 4.DO ENDKINCO

BEGIN EXTERNAL LOAD 50.DO

END EXTERNAL LOAD

BEGIN STIFFNESS CONSTANTS PER LINE ELEMENT

1

.D0,2.D0 9.DO12.DO

1 .DO,%.DO