“k'é'fjit—e ' -

. . n
wan. .
2",: ‘3»
.‘w ’~ _

.\.

 

' A.» L. t4 “qu1-

I.’IZ'. M.
r: V“:
”5;“
.. ‘

KNEE.

Mi ‘ ”a « 192'."
, u.;.-\\ 3 . . 9L‘.;.‘.
‘- , l» .1! x ‘ ' ’ a '
_ .‘.‘..,.‘..'.¢-.’~.t ‘ "1 ‘ y. .51.:
‘ ' u 5&1»; (
.
3?; V . ~ "
" . “1‘1“? 'L 1:- ' -
. ' ‘ ‘1» , u .2. ww-
"!1‘§5“3u.~ .m . 7‘ 7 many: -= 1." k
’v‘fifi j*;v.‘»‘3.'ww‘3~ “r ’1‘" ,— gram. .-
.-- ' ““' V 1‘ ..n w“
--u‘~-*‘-'r;.°:z:a a r {"3}! i ». “.0" ”Exam:-
, n ‘ . 134- u
» :r

*2; mm
. >> M say) a??? J» .
an! vi -‘ ' ”if. “fir. ”L .w..r»1
bflrAn-a ‘ 3'. . orig“ - I‘nr
p

. ,3:

. "t
..~:g:ȣ.m
HM. .. .
‘1’"1‘5':

‘ b
M" ."I A.
r: z ‘ v r: 5r.

fo-v1-YT

_. . t ‘ 'ht’kw
c. . - V . . . l . ,. ‘ ~
- 3- . ‘ . a. ‘ MW” . , ‘
. , n “MAMA“. ..x . -,. Tint?- , ‘
-‘ HR-v."“«.".‘-."" n A ,. (V:- (g ,. ”.55 'v ., 7 V , ,
5.. ‘35 1t 3, .., V , ,3" "3 5” “7:75.,“ ‘ ‘ . “H > . 1‘“;
\.- ‘ ‘ A‘JW‘.‘ ‘ ":3? ti- ’- “A $ '1 ., _ -' w
’ J" ‘- " l“ I}! X: A 1' ' W i ', ' t - ' ‘ 4 1.35:. it
. L} V’ ‘1 I . ._ 0‘» .‘E 5 v 'I 5'33,“ Ff. ‘ ’EL . H , f'n , ref", Sfim-
f ‘ Y .3.- ' , V. . , > _ ’ .. "Jr:
,r ‘ h. ' ”A - «r;- - .
A . ”#4
‘ A

~c;-: .
h.) J" “'
:44 kg»;

.|04—

 

.m‘w

9.:-
.

b' ‘v
.95.? _
‘3.\“ n1. *3 ,
fizzy.» 4., q‘
w‘hn M
a): .1» n ;,
w grub». w." 3 ~ , 7 'u ,-
4!» I , ‘ A > ' .
" ' arms" 95:? 3'
"In": ugly-ca . mam“?
4‘ «arm: 1

 

THY—:51?

l inn Will/Iii » i =~ ~

079476 P’ '5

“N
Eat
\'\\
.. 7:"
t " “\‘u
.. ‘ 9‘s.»
\ ‘
val-nu"-

? ' tier: 2:: say

«m! I

This is to certify that the

thesis entitled

Modelling Cylindrical Discontinuities
In Running Shoe Soles Using A
Condensed Finite Element Formulation

presented by

Andrew James Hull

has been accepted towards fulfillment
of the requirements for

 

 

   
  

    

M.S. degnfigh] Mechanical Engineering
Major professo’ ' '
Clark Radcliffe
[hue February 21, 1985

 

MS U is an Affirmative Action/Equal Opportunity Institution

 

 

MSU RETURNING MATERIALS:
Place in book drop to
LngAxlgs remove this checkout from
w your Y‘ECOY‘d. FINES_ will

be charged if book is
returned after the date
stamped below.

 

 

 

 

 

Q {@214 U 4
“ ‘ 4%;e
W

NOV 16 2004

 

 

MDELLING CYLINDRICAL DISCONTINUITIES 1N RUNNING SHOE sows
USING A CONDENSED FINITE ELEMENT FORMULATION

By

Andrew James 3111 l

AmESIS‘

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Mechanical Engineering

1985

3346/33

ABSTRACT

MODELLING CYLINDRICAL DISCONTINUITIES IN RUNNING SHOE SOLES

USING A mNDENSEn FIN ITE ELEMENT FORMULATION
By

Andrew James Hull

This thesis presents modelling techniques that are useful to
analyze the effects of holes in running shoe soles. A.two dimensional
finite element analysis of a running shoe containing one to three
horizontal holes in the heel area was performed using condensed two
dimensional elements. The displacement results serve as a basis for
the effects of these holes during heel strike. A three dimensional
element program that condenses a twenty node solid into a ten node
solid and formulates the element stiffness matrix was written. The
practical and economical advantages of using condensed formulations

were demonstrated.

DEDICATION

To my father

ii

ACKNOWLEDGEMENTS

I would like to express my sincere appreciation to my major
professor. Dr. Clark Radcliffe. His knowledge. guidance, and
friendship have made this research possible.

I would like to thank.Wolverine World Wide. Inc.. for funding
this research and Swanson Analysis Systems. Inc.. for the use of
ANSYS.

' Special thanks to my mother. brother. and sister for their
constant support they have given me the entire time I have been in
school.

I would also like to thank Drs. R.W. Sotus-Little. M.v.
Gandhi. and Erik Goodman whose comments and advice were greatly
appreciated.

Finally I would like to thank Ken Rowe. Ken Specht. Brian
Sterner. and Dan Larabell for their friendship during my academic

career 0

iii

LIST OF TABLES . . . .
LIST OF FIGURES . . . .
NOMENCLATURE . . . . .
INTRODUCTION . . . . .
TEE FTNITE ELEMENT'METEOD
TWO DIMENSIONAL ANALYSIS .
THREE DIMENSIONAL ANALYSIS
CONCLUSIONS' .' . . . .
REFERENCES . . . . .
APPENDIX A . . . . .
APPENDIX B . . . . .
APPENDIX C . . . . .
APPENDIX D . . . . .
APPENDIX E . . . . .
APPENDIX F . . . . .

APPENDIX G . . . . .

TABLE OF CONTENTS

iv

Page

. 20

. 27

. 30

. 46

O 75

Table

Table

Table

Table

Table

Table

Table

Table

Table

Thble

Table

2.1
2.2
3.1
3.2
3.3
F.l
F.2
F.3
F.4
F.5

F.6

LIST OF TABLES

 

X displacements for Figures 2.1. 2.2. and 2.3 . .

Y displacements for Figures 2.1, 2.2. and 2.3 . .

Material pr0perties of the Brooks Supervillinova .

Displacements in the x direction

Displacements in the y direction

Y direction displacements
Y direction displacements
Y direction displacements
X direction displacements
X direction displacements

X direction displacements

of

model

model

model

model

model

model

with

with

with

with

with

with

one hole .

two holes .

three holes

one hole

two holes .

three holes

11

15

15

77

77

78

78

78

79

Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure

Figure

Figure

Figure
Figure
Figure
Figure
Figure

Figure

Figure

2.1
2.2
2.3
3.1
3.2
3.3
3.4
3.5
5..
3.7

3.8

3.10
3.11
A.1
3.1
8.2
3.3

3.4

Figure 3.5

Figure

3.6

LIST OF FIGURES

Base grid with condensation . . . . . . . .
Base grid without condensation . . . . . . .
Refined grid with linear elements . . . . . .
Brooks Supervillinova running shoe . . . . .
Base grid of a Brooks Supervillinova . . . . .
Average vertical component of ground reaction force
Average shear component of ground reaction force .
Refined grid of a Brooks Supervillinova . . . .
Displacements of node 5 versus radius size . . .
Displacements of node 5 versus radius size squared
Displacement plot of sole with one hole. R!.30 .
Displacement plot of sole with two holes. RP.30 .
Displacement plot of sole with three holes. R!.30
Out of plane shear force . . . . . . . . .
Tho dimensional triangular element . . . . .
Four node two dimensional solid element . . . .
Shape function of a four node element . . . .
Eight node two dimensional solid element . . .
Shape function for an eight node element . . .
Eight node three dimensional solid element . '. .

Twenty node three dimensional solid element . .

vi

10
11
12
13
14
16
17
18
18
18
19
36
47
49
49
50
51

52

Figure 8.7 Five node two dimensional solid element . . . . . 54

Figure 3.8 Four node element mapped into natural coordinates . 58

vii

NOMENCLATURE

mapping matrix. x direction
mapping matrix. y direction
mapping matrix. 2 direction
constraint matrix

stiffness matrix

shape function

modified shape function
natural coordinate position
natural coordinate position
natural coordinate position
generalized position
spatial coordinate position
spatial coordinate position
spatial coordinate position
generalized displacement
natural coordinate r. s. or t
stress component

shear stress

viii

INTRODUCTION

Shoe designers frequently place holes in running shoe soles. The
orientation of these holes is both vertical and horizontal in the shoe
sole. The number of holes also ranges considerably: some shoes do
not have them. and other shoes incorporate more than 35 holes into the
sole design. There are various reasons for including these holes.
ranging from weight reduction to aiding in manufacture of the shoe.
The effects of these holes has never been fully understood.

This thesis applies a finite element analysis to a two
dimensional shoe sole containing one to three holes orientated
horizontally in the wedge area. The radius of these holes is also
varied to study the effects of different size holes. Due to the
number of nodes required, the current finite element pragram in use
does not have the ability to model the three dimensional problem. A
three dimensional modelling technique is discussed using nodal
condensation to reduce the three dimensional problem to a manageable

size.

TEE FINITE ELEMENT’METEOD

Analytical methods in structural analysis have been studied for
many years. Exact solutions are usually not available to problems
with complicated geometry and/or boundary conditions. The emergence
of high speed computers. however. makes it possible to accurately
apply numerical methods to complex problems. A commonly used
numerical technique for structural analysis is the finite element
methodl1.8]. This method involves a discretization of a continuous
structure into a number of smaller parts (elements). Equations are
first formulated for each element with individual loading and boundary
conditions. then a set of equations are assembled modelling the entire
system. The resulting equations are then solved. yielding the
structural response.

For an elastic body. the stiffness matrix and load vector can be
formulated using the finite element method. and the displacements of
the body can be determined (Appendices A. B, C, and D). The best
'method to verify the results of any analytical technique is actual
testing. When testing is not feasible. a commonly used convergence
test is to reduce the element size and compare the results of a denser
grid of smaller elements to the base model. If the differences in
responses are small. then the base model probably has enough elements
to model the system. This was the verification technique used in this

paper. comparing solutions of variously populated grids to insure

analytical convergence.

The major tapic of this thesis is the effect of holes in running
shoe soles. therefore. a modelling technique was develOped to
determine the global effects of holes in materials. Consider the
three 2 dimensional element test grids shown in Figures 2.1. 2.2. and
2.3. All three were subjected to a single point load on the upper
right node. The four left hand nodes were constrained in the x and y
direction. The structure was six units long. six units high. and the
hole in the middle was one unit in diameter.

There was a choice of two dimensional models based on modelling
assumptions. The model assumption of plane strain was used. Plane
strain is a specialization of three dimensional linear elastic theory.
It represents the situation where the structures geometry and loading
are constant in the z direction and the component of displacement
perpendicular to the r-y.plane is zero. This is discussed in Appendix
E. Plane stress. the other two dimensional specialization of three
dimensional elastic theory. represents the situation of a very thin.
flat structure. whose loading occurs only in the x-y plane of the
plate. All three models were deve10ped and analyzed using ANSYS.
ANSYS is a general purpose finite element program developed by Swanson
Analysis Systems. Inc.

Figures 2.1 and 2.2 illustrate the five and eight node elements
used to model the hole in the middle of the soiid. Each element
bordering the hole had three active nodes on the hole boundary. The
midside node on the boundary of the hole was left on to insure a curve

was fitted thru this area for the element interpolation. Different

13

 

 

 

 

 

 

 

 

 

 

 

 

>

 

'X

Figure 2.1 Base grid with condensation

 

 

 

>

X

Figure 2.3 Refined grid with linear elements

sources recommend different locations for this node[1,s], From a
theoretical standpoint. the midside node must be located near the
middle of the curve between the two end points. When the midside node
is too close to the end nodes of the curved side. the Jacobian matrix
becomes singular. and has no inverse. The inverse of the Jacobian
matrix must exist to insure a unique mapping from spatial coordinates
to natural coordinates (Appendix B).

Figure 2.1 was the base model. It consisted of 36 elements and
32 nodes. The four elements in the middle of the model were five node
isOparametric solid elements. They each had three nodes condensed out
of the element matrices (Appendix B). Figure 2.2 was the second base
model. consisting of 40 elements and 40 nodes. The four elements in
the middle of the model were eight node isoparametric solid elements.
This model was included so that the effects of nodal condensation
could be viewed. Figure 2.3 was the refined model. It consisted of
108 elements and 100 nodes. All the elements of Figure 2.3 used
linear interpolation. The results were compared in Tables 2.1 and
2.2. The displacements of the base models in Figures 2.1 and 2.2 were
compared to the displacements of the refined model in Figure 2.3.
Nodes l thru.8 were on the circle. and nodes 9 thru 12 were labeled on
each figure. Node 13 and 14 were the right corner nodes. node 13
being located directly under the load.

The results from.Thble 2.1 shows that the local convergence (the
displacements on the hole) in the x direction was very poor. The
grids in Figures 2.1 and 2.2 would not be adequate if the

displacements on the circle were desired. The global convergence (the

 

 

Table 2.1 X displacements for Figures 2.1, 2.2. and 2.3

Node Figure 2.3 Figure 2.1 Normalized Figure 2.2 Normalized
number -refined grid -base grid difference -base grid difference
a b b

1 ’.0029 -.0931 -2.5 -.0698 -1.9
2 -.0651 -.1863 -3.4 -.1742 -3.1
3 -.1735 -.2198 -1.3 -.2306 -1.6
4 -.2225 -.1816 1.4 -.2101 0.3
5 -.2172 -.1183 2.8 -.1485 1.9
6 -.1961 -.0682 3.6 -.0887 3.0
7 -.1464 -.0713 2.1 -.0940 1.5
8 -.0533 -.0636 -0.3 -.0627 -0.3
9 ‘.2631 -.3071 -1.2 -.2732 -0.3
10 -.6614 -.6891 -0.8 -.6689 -0.2
11 .2892 .3262 1.0 .2997 0.3
12 .1864 .2239 1.1 .1990 0.4
13 3.5680 3.5370 -0.9 3.5590 -0.3
14 ~2.2870 -2.2630 0.7 -2.2810 0.2

Table 2.2 Y displacements for Figures 2.1. 2.2. and 2.3

Node Figure 2.3 Figure 2.1 Normalized Figure 2.2 Normalized
number -refined grid -base grid difference -base grid difference
a b b

1 -1.858 -1.890 0.4 -1.889 0.4
2 -2.486 -2.429 -0.7 -2.455 -0.4
3 '3.104 -2.951 -1.8 -3.005 '1.2
4 -3.378 -3.174 -2.4 -3.232 -1.7
5 -3.112 -2.990 *1.4 -3.018 -1.1
6 -2.455 ~2.450 -0.1 -2.431 -0.3
7 -1.822 ~1.896 0.9 -1.855 0.4
8 -l.584 -1.682 1.1 -l.673 1.0
9 -l.283 -1.288 0.1 -1.275 -0.1
10 -3.758 -3.635 -1.4 -3.712 -0.5
11 -3.885 -3.745 1.6 -3.842 -0.5
12 -l.200 —1.196 ~0.0 -1.196 ~0.0
13 '8.590 -8.424 -1.9 -8.532 ~0.7
14 -5.578 -5.433 -1.7 “5.525 -0.6

Normalized difference = [(a1 - bi)/max(a)] X 100*

displacements at common points not on the hole) in the x direction was
fair in Figure 2.1 and good in Figure 2.2. The results from Table 2.2
shows that the local and global convergence were good to excellent for
Figures 2.1 and 2.2. For the two dimensional model. the modelling
technique in Figure 2.1 was used. For each hole modelled the number
of nodes needed using a five node element was eight. For the
conventional method of using many linear elements. the number of nodes
needed was 68. Since each node had two degrees of freedom (adding two
equations per node to the system of equations) the difference in
number of equations between the base grid and the refined grid was
128. Furthermore. the formulation time for four 5 node elements was
less than the formulation time for 64 four node elements. The use of
higher order elements in this case resulted in a time savings on the
computer while retaining the numerical accuracy in the the finite

element analysis.

TWO DIMENSIONAL ANALYSIS

The first analysis of the running shoe was a model of a Brooks
Supervillinova (Figure 3.1). A two dimensional model of the
heel/wedge area was developed (Figure 3.2). which consisted of 82
elements and 103 nodes subjected to a plane strain analysis. The shoe
in this model was cut with a plane that ran parallel to the longest

direction of the shoe and was perpendicular with the ground.

 

Figure 3.1 Brooks Supervillinova running shoe

The bottom row of elements in Figure 3.2 represented the outsole. the
middle three rows of elements represented the midsole, and the top row

of elements represented the sockliner. The sockliner was a removable

10

11

layer of material placed in the running shoe to increase comfort. The
model was constrained in the x and y directions at the nodes marked
with triangles. Although there was a tread pattern on the outsole.
this area was modelled as a solid since the protrusions occupied over
85% of the volume under the loaded area. The material preperties were
determined using a Model 1331 Instrom Servohydaulic Materials Testing
Machine located in the Biomechanics Department at Michigan State

University. They are listed below in Table 3.1.

Table 3.1 Material properties of the Brooks Supervillinova

 

 

 

Layer Young's Modulus Poisson's Ratio

Outsole 260 psi .25

Midsole 130 psi .25

Sockliner 80 psi .25
“m
\K:\‘\.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

L.

 

 

Figure 3.2 Base grid of a Brooks Supervillinova

The section of the size eight and a half shoe modelled (Figure 3.2)

was 5.2 inches long and 1.2 inches thick. Shoe designers who include

 

12

horizontal holes in their midsole design normally locate them in the
heel area of the shoe. The model was developed to investigate the
effects of these holes.

The following assumptions were used: 1) All sole materials
behaved elastically. 2) the shoe upper provided negligible additional
stiffness to the sole 3) there was no slip between the ground and the
bottom of the shoe. 4) The runner struck the ground with the heel area
of his shoe. 5) Plane strain occured.

The loading for this model was derived from experimental

measurements of pressures exerted by running human subjects.

 

 

 

50 10b 150 2330 TIME

Figure 3.3 Average vertical component of ground reaction force

iFigure 3.3[4] was the vertical component of ground reaction force in
1L2 subjects. running at 6.5 minutes per mile. The solid line on the
Staph was the average value and the dashed lines were the maximum and
'minimum ranges. Units of force (Fr) were in body weight and the units
03 time were in milliseconds. The ground reaction force in the x

dii-xection was vertical. and a positive force was orientated in an

13

upward direction. The first maxima of the curve had a value of 2.1
'units of body weight at approximately 25 milliseconds. Since this
finite element model only analyzed wedge area displacements. the first
local maxima was used as it corresponded to heel strike. Using a
subject weight of 150 lbs. the reaction force at heel strike in the
shoe was (-2.1 X 150)= -315 lbs. Since the width of the shoe was 3
inches. and a plane strain analysis has a unity thickness. the loading
for this model was (-315 / 3)= ~105 lbs total in the x direction.
Distributing the -105 pound lead over ten nodes resulted in a force
vector at each node equal to (-105 / 10)= -10.5 lbs per node in the x

direction.

 

Y // \\
;//\ \\
P \ 31* 1
1 \ F ' 7 '
\ \\\// //T|ME
/
\\’//

 

-1 4—

Figure 3.4 Average shear component of ground reaction force

Figure 3.4[‘] was the average shear component of the ground
reaction force. The ground reaction force (Fy) in the longitudinal y
direction was determined using Figure 3.4. Positive force was
f01rward. At 25 milliseconds the average shear component was .2 units
°f bOdy weight. Using the same method of calculation as above. the

everage shear force per node was equal to 1.0 lbs in the y direction.

14

The ten loaded nodes were denoted with vectors in Figure 3.2. For
convenience they were called nodes one thru 10. one being the left
most node. The analysis was also run using a finer grid of 326
elements and 346 nodes. shown in Figure 3.5. The results are listed
in Tables 3.2 and 3.3. With the exception of the left most node. the
y direction results were fairly close. The y direction displacements.
or the amount of compression of the sole. is the critical response in
the design of a shoe. This indicates that the base grid (Figure 3.2)
was fine enough to model the load and constraints that were being
applied to it. A grid modelling the entire sole was also analyzed.
using the constraints and loading conditions determined above. The
percent difference for the displacement of all common nodes was less

than one percent. This justified ignoring the forward part of the

shoe as in Figures 3.2 and 3.5.

  

L.

X

Figure 3.5 Refined grid of a Brooks Supervillinova

‘Using five node two dimensional elements. the grid in Figure 3.2
Was modified to include the effects of holes running across the

midsole from the medial to the lateral part of the shoe. The

Node
number

o‘omucxM-waD-l

H

Node
number

c.m>o>~ao~u~e.oaesr¢

H .

Displacement
'-base grid

b

.01564
.00482
.01540
.02995
.04409
.05633
.06533
.07061
.07658
.05822

Displacement
'-base grid

b

“.2992
“.2909
“.2923
“.2867
“.2836
“.2819
“.2807
“.2779
“.2691
“.2427

15

Normalized difference

Displacement
-refined grid

.02456
.02050
.02712
.03743
.04843
.05842
.06637
.07152
.07643
.06437

Displacement
-refined grid

“.2268
“.2754
“.2836
“.2840
“.2826
“.2815
“.2804
“.2778
“.2716
“.2364

Table 3.2 Displacements in the x direction

Normalized
difference

P‘bfih‘
\Otncssa
O

I
OOOHHNu
c>ezhae-l~l¢3usUIa\

Thble 3.3 Displacments in the y direction

Normalized
difference

I
N

éééoonbua
Laboésleiainioioioeo

i
N O

[(ai — bi)/max(a)] X 1005

16

displacement effects of one. two. and three holes were determined
using radius sizes ranging from .05 inches to .30 inches. incremented
by .05 inches. All of the holes had their center located at y=.55
inches from the origin shown in Figure 3.2. For the single hole. the
location was at x=2.05 inches. For the analysis using two holes. the
x locations were at 1.45 and 2.35 inches. The model that included
three holes used x locations of 1.15. 2.05. and 2.95 inches. A
displacement plot of a typical node is shown in Figure 3.6. The
tabular data for the nodal displacements of nodes one thru 10 for all

geometry modifications is listed in Appendix F.

UY
-060 '-
————— THREE HOLES

-.55- ---------- TMOIKAES }D
-.58
-.45
-.48

‘035

 

‘03“.

 

l l l l l L

.85 .18 .15 .28 .25 .38 R

Figure 3.6 Displacements of node 5 versus radius size

.A general trend of the data indicated that as hole radius or the

number of'holes were increased. the absolute value of the displacement

 

17

increased (Figure 3.6). For the geometry range tested. the
displacement was roughly proportional to the area of the hole (Figure
3.7). The nodal displacements were also sensitive to the location of
the holes. Nodes above the holes deformed to a greater extent than

the nodes above solid material (Figures 3.8. 3.9. and 3.10).

UY
' 060 '-

——————— THREE HOLES
- .55 - ------------- Two HOLES /
ONE HOLE // ,A

 

"050

-045

-040

.035

-030

 

 

L l l l l l

. 8188 . 8225 . 8488 . 8625 . 8988 RNZ

 

Figure 3.7 Displacement of node 5 versus radius size squared

The main advantage of the two dimensional analysis was rapid
solution time. when compared to the three dimensional model. because
the grid was simpler and] fewer nodes were present in the model. This
modelling technique however. was a two dimensional approximation to a
three dimensional problem which had finite depth. hence plane strain

did 'not occur. Furthermore. some shoes have the holes running in the

18

 

 

 

Figure 3.8 Displacement plot of sole with one hole. R=.30

 

 

 

 

Figure 3.9 Displacement plot of sole with two holes. R=.30

 

 

 

 

Figure 3.10 Displacement plot of sole with three holes. R=.30

 

 

 

vertical direction which cannot be modelled using a two dimensional
approach. The two dimensional approximation agreed with measured
values because the out of plane shear force in the z direction was

very small and varied widely from one runner to another (Figure

 

 

 

3.11141),
0.5 "I
F2 /<;:;>r_.‘~/’ ”:-_"—'=é§;553===:=,. L
\_’ ’ \;__\\:’:’: / _ ' TIME
-O.5 "

Figure 3.11 Out of plane shear force

THREE DIMENSIONAL ANALYSIS

The three dimensional analysis of structures is more complex than
the two dimensional counterpart. The extra degree-of-freedom at each
node causes the material stiffness matrix to enlarge from 3X3 to 6X6.
The stress components in the finite element formulation changes from

161T - fox oy txy} (4.1a)

to

{air = {orx o

y o t t t l . (4.1b)

The element considerations in the finite element fonmulation also
changes. Because three dimensional curved side elements require
midside nodes like their two dimensional counterparts. the definition
of a three dimensional curved side element with midside nodes requires
20 nodes instead of the eight present in the two dimensional curved
side element. The wave front limitation of the ANSYS educational
version is 200 and the connectivity of 20 node three dimensional
elements prevent an analysis of a three dimensional model that
incorporates vertical holes.

Condensation of the current 20 node solid into ,a smaller (in

terms of number of nodes) element is a technique that can reduce this

20

21

problem to a size where an analysis is possible. Quadratic behavior
on one side of the condensed element is desirable so the element face
can follow a curved surface such as a hole. Curved boundaries cannot
be followed by linear elements. The eight midside nodes that do not
lie on the hole can be removed through condensation because they do
not lie an curved boundaries. The additional nodes that lie on the
hole can be condensed because the holes under consideration are
constant radius and do not have a curve in their axial direction. The
result is a ten node condensed element that can incorporate a quarter
of a hole on one side as a boundary condition. The quadratic
displacement behavior along the edges where condensation occurs
becomes linear. This is not undesirable since there is almost no loss
in accuracy.

Linear three dimensional elements attach to the condensed element
in a much more efficient pattern. This is due to element mesh
compatibility of the condensed sides of the ten node element with any
side of the eight node three dimensional linear element.

The condensation procedure started with a twenty node finite
element map to establish the Jacobian matrix. The 20 node element
equations that defined the mapping from natural coordinates (r.s.t) to
spatial coordinates (X.Y.Z) were '

X = a1 + agr + a,s + "t + a,r2
+ tgsz + t7t2 + e.st + a,sr + a,.rt

2

+ ‘1132t + 913323 I ‘ast23 + ‘aatzr I ‘asr 3

+ a1‘:2t + a1,rst + a1.32rt + a,,t2rs + a3.r28t (4.1a)

 

22

Y =b1 + bzr+ b38+ b4t+ bsrz

+ b :2 + b t2 + b.st + b sr + b .rt
‘ 7 9 1

+ b11s2t + blzszr + b13t2s + bljtzr + b15r2s

+ b1‘r2t + b17rst + b1gs2rt + b1,t2rs + baorzst (4.1b)

- 2
Z - c1 + c,r + c,s + cgt + c,r

2 + c,st + c,sr + exort

2 2

+ c.s2 + c,t

2

+ caa‘zt + “123 r + °sst 3 + c1‘t2r + °isr s

2rt + c1,t2rs + c,,:2st (4.1s)

+ casrzt + °1133t + cis‘
The vectors a. b. and c defined the mapping and were constant for each
element defined. They were found by mapping the element back into
natural coordinates. which produced three systems of twenty linear
equations. the solution being the a. b. and c vectors. This process
is similiar to the two dimensional mapping discussed in Appendix.B.

Tb remove midside nodes and condense the element. the shape
functions were modified in their natural coordinates. This step
differs from most finite element formulations. The shape functions

defined by the 20 node element were condensed using[91

9m = (1/2)(“c1 + “c2) (4.2a)

where n represents r. s. or t. the subscript m denotes a condensed
midside node. and the subscripts c1 and c2 denote the corner nodes

that are collinear with node m. Applying this equation at the ten

midside

nodes

direction from

20
8(r.s.t)=§1Nj(r,s.t)8j

10 .
5(1’.S.t) = N'( a Dt)6'

where

5

is either r.s.

N1 +
NJ +
Nx +
NL +
NH +
NN +
N0 +
NP +
No

NB 0

(NT/2)
(NR/2)
(NR/25
(NS/2)
tux/2)
(NV/2)
(NV/2)
(NW/2)

+

+

+

+

+

+

+

+

changed the

or t.

(NY/2)
(NZ/2)
(NS/2)
(NT/2)
(NY/2)
(NZ/2)
(“w/2)
(Nx/z)

23

displacement field approximation in the 6

(4.2h)

. (4.2c)

The modified shape functions became
(4.2d)
(4.2c)
+ (NA/2) (4.2i)
+ (NB/2) (4.2g)
(4.2h)
(4.2i)
+ (NA/2) (4.2j)
+ (NB/2) (4.2k)
(4.21)
(4.2m)

This modification eliminated ten nodes per element. Each node had

three degrees

of freedom.

therefore. the condensation technique

eliminated thirty equations per element. This technique is discussed

in Appendix B and an example of a two dimensional element condensed at

24

three locations is included.

The derivatives of the shape functions were next multiplied by
the inverse of the Jacobian matrix. This step is described in
Appendix B. The Jacobian matrix and it's inverse were defined using a
twenty node map. while the derivatives of the shape functions were
defined using a ten node condensed element. The stiffness matrix was
then defined using the equations from Appendix B and numerically
integrated using a 3X3X3 grid and the Gauss-Legendre quadrature method
(Appendix C).

The element program in Appendix G will build the element
stiffness matrix for ten node degraded elements. To insure pregram
accuracy. the stiffness matrix of an ANSYS defined twenty node element
was output and this stiffness matrix was modified to incorporate the
linear constraints implied by the condensation process. A constraint

matrix was defined by
{u} = [Cliul' (4.3a)

where u is the vector containing the x. y. and z displacements of the
60160 element matrix. u' is the x. y. and z displacements of the 30X30
matrix. and C is the' constraint 'matrix (60X60) that contains the
linear displacement constraints on the midside nodes. The constraints
were applied to the ANSYS stiffness matrix yielding a new system of

equations. as shown in equation (4.3b).

mm = [eiTmtcuui' =- [Kl'Iul' . (4.31.)

 

25

This is best illustrated by the three by three system of equations
shown below
k11 k1: R1: “1
kn kn k3: n. = [xltul (4.4a)
kll ks: kss “a
which is constrained such that
u, = n; (4.4c)
u, = (mm; + u;) . (4.4.1)
The constraint matrix becomes
1 0 0
[C] = .5 0 .5 (4.4a)
0 0 1
and assembling this according to equation (4.3b) produces a modified
stiffness matrix
k11+a5k31+e5k13+e25k13 0 k13+05k13+05k23+'25k33
[X] 8 0 0 0 (4.4f)

k31+.5k3;+.5k33+.25k33 0 k33+.5k31+'5k33+°25k33

26

Using this method. the results of the condensed element pregram were
compared to a stiffness matrix multiplied by constraint matrices. The
percent difference averaged 0.11% at all nonzero locations in the
matrix. This method has the advantage of producing a stiffness matrix
to compare to the ten node condensed element stiffness matrix without
running an actual finite element analysis.

There is considerable advantages using condensed elements in
three dimensional analysis. A.running shoe sole such as the Brooks
Supervillinova has 37 holes running vertically. The number of
equations for this problem without condensation is approximately
36000. which is too large for almost any finite element routine.
Since three layers of four elements are needed to model each hole. 544
condensed elements are needed to model this sole. Each condensed
element eliminates thirty equations. and the total savings for this
type of analysis is 31320 fewer equations. The size of the new
problem using condensed elements is approximately 4680 equations. or
roughly a magnitude smaller problem. The disadvantage is that the
quadratic behavior along the boundaries becomes linear. This.
however. is usually not a concern as long as the finite element grid
is distributed in a manner that several linear elements fit the
displacement curve. Unfortunately ANSYS does not currently have the
ability to accept elements input by the user and no actual analysis

was performed using this element.

CONCLUSIONS

This thesis has shown that using a condensed finite element
formulation can greatly reduce the problem size of three dimensional
models. The numerical accuracy using two dimensional condensed
formulations was retained. as demonstated by a two dimensional
analysis. A program that generates ten node condensed element
stiffness matrices is available now to be interfaced with finite
element programs that can accept user defined elements. The accuracy
of this element was compared to an ANSYS constrained 20 node element.
and the percent difference was found to average 0.1rs at all nonzero
locations. The element has the advantage of eliminating thirty
equations every time it is used instead of the base three dimensional
element.

The displacement effects of horizontal holes in the sole of a
shoe were calculated. The analysis included one to three holes
ranging in radius from .05 to .30 inches in the wedge area of the
shoe. No attempt was made to determine if these displacements were
favorable or detrimental to the runner.

Factors not addressed in this thesis are the dynamic effects of
the shoe material. the dynamic effects of the loading. the effect of a
ground strike in a location other than the heel area. and the
interaction of the shoe upper with the sole. Future tepics also

include development of a program that uses the three dimensional

27

28

condensed element stiffness matrices. assembles them along with the
load vectors and linear three dimensional elements. and solves for the
displacements of the nodes. resulting in a finite element routine to

handle large scale problems.

REFEREI CES

REFERENCES

 

Zienkiewicz. O.C.. Th3 Finite Element Method, McGraw-Hill Book
Company. 1983. '

Ruebner. Kenneth H.. Ihp Finite Element Mpthod jg; Engineerg.
John Wiley and Sons. 1975.

 

Kohnke. Peter C.. ANSYS Engineering Analysig System Theorptipgl
Manual. Swanson Analyis Systems. Inc.. 1983.

 

Cavanagh. Peter R.. "The Shoe-Ground Interface in Running".

Americgn Apggpgp pf Orthopaedic Surgpons Sppppgium pp Th; Fppp
pp; Lpg ip Running Sports. The C.V. Mosby' Company. 1982. pp.
30-44.

Cook. Robert D.. Qoncepts ppg Application; pf Fipipe Elemept
Apalpsis. John Wiley and Sons. 1974. '

Carnahan. B.. Luther. H.A.. and Wilkes. 1.0. Applipd Numeripal
Methods. John Wiley and Sons. 1969.

29

APPENDICES

 

APPENDIX A

FORMULATION 0F ms STIFFNESS MATRIX AND LOAD VECTOR

FROM THE MINIMUM COMPLEMENTARY ENERGY PRINCIPLE

The stiffness matrix corresponding to an elastic body can be
derived in a number of different ways. One way to formulate this
matrix is by starting with the complementary energy principle. which
corresponds to a compatibility condition. The principle states that
the state of stress that satisfies the stress-strain relations in the

interior and all prescribed boundary conditions also minimizes the

system's complementary energy in an elastic body. If
“c1“x'°y'°z'txy'txz"yz? is the complementary energy.
Uc(ax'ay'°z'txy’txz'tyz) is the complementary stress energy. and Vc is

the work done by the applied loads during stress changes. then.

according to the minimum complementary energy principle.

bnc = 6(Uc-Vc) (A.1)

= 60° - avc = o . (A.2)

Tb facilitate a finite element formulation. the variation is taken

30

31

with respect to the stress components. The complementary stress

energy is defined as

Uc(ox.ey.ez.txy.txz.tyz) = l/ZJIIvfoITIDlioldV (A.3a)

where

y 2 xy txz tyz} (A.3b)

161T 3 {ex 6 o t
in a three dimensional problem. The matrix [D] is the material
stiffness matrix that is used to relate the strains to the stresses.
In a three dimensional problem it is a 6X6 matrix. In a two
dimensional problem it is a 313 matrix and has different values
depending on if the problem definition is plane strain or plane

stress. It is defined by
{c} = [D1161 (A.3c)
where the strain vector

MT = {a (11.3.1)

x ‘y 8z 7xy 7xz Vyz} '

The relationship between the stresses and the strains in equation

(A.3c) can also be written as

{o} = {Clio} (A.3e)

 

32
where
[c] = [01'1 . (A.3f)

Including the column vector of initial strains is.) into the

complementary stress energy term results in

n. = 1(2m.,[taiTin1iai — 2ia}T{..i]av . (11.4)
The work done by the external forces is

Vc . ”Lb.“ + r‘v. + z‘w.]dv + fl,[r;u. + 1;“ + T;w.]dS (A.5a)
or written in matrix terms

vc .. jflvuzfi'rrsidv + ILH‘libldS (11.51.)
where {FIIT are the body forces such as gravity. [TI] are the boundary
traction components acting on the surface of the solid. and {6} is the

column matrix of the components of the displacement field.

Substituting equations (A.5b) and (A.4) into equation (A.2) results in
1’2111v1“’T[”1{‘i’Z‘“’T“°’1“”.111v‘F-‘fl‘5mv'fls“.1{WSW (A. 6)

Substituting equation (A.3e) into (A.6) yields

33

1/2Ifjv[(rc1(.iiTinitciia} - 2[[Cl(s}lT{c.}]dv
- ”Lue‘fl'isidv - ”,[r‘ltsms = ac . (Ac?)

Now defining the relationship between the strains and the

displacements as
{a} - [31(6) (A.8)

where [B] is a 3NX6 matrix (for a three dimensional solid with N nodes
per element) that contains the derivatives of the shape functions.
Appendix B discusses the shape functions and their derivatives. Using

equation (A.8) in equation (A.7) produces

1(2IIIV[[[31(611TICJTIDJ[cits]is} - zltslieiithiTie.1]av
- ”Ltp‘itsidv - ”'[Ifitsias - no . (A.9)

Since the material properties are isotropic. the material stiffness

matrix is symmetric. i.e..
[c] . [CIT . (A.10)

Substituting equations (A.3f) and (A.10) into equation (A.9) will give

the following result

ilzfjjv[ts}TIBIthilsiis} - 2(siTrsiTrcir.,}]dv
- j'ILIF‘iTtsidv - ”.[r‘ltsids - no . " (5,11)

34

For a single element (e) the discretized functional is

fl(°)({8}(°)) 3 “(e)(u1.u3.....ur.v1.vz.....Vr.w1.'3.....wr) (A012)
or in terms of equation (A.11) the functional is
“(e)({5}(e)) g 1/2IIIF[{8}T(°)[B]T(°)[C](°)'5}(°)
- 2{5)T(e)[B]T(e)[c](e)(8.}(e)]dv(e)
.. IIIV{F.}T(°)'5} (O)dv(e) .. II8[T’](O){5}(O)ds(¢) . (A.13)

At equilibrium the potential energy of a system is stationary. The
discretized system is stationary when the first variation of the

discretized functional vanishes. thus
I (e)
5n(u.v.w) 8 2 [6n (u.v.w)] 8 0 (A.14a)
. . 931
where

6n(’)(u.v.w) 8 2:'1[(3n(°)laui)6ui]

-1 . . -1 .
(e) ‘ (e)

+ 2i-1l(an lavi)6vi] + 21'1“” la'iis'i] , ”.141”

Th0 531. bvi. and the bwi are independent (they may or may not be

zero) therefore. the individial parts of the summation are forced to

zero. i.e..

«kW/ani = “(Q/avi = 3n(°)/8wi = o i=1.2.....r (A.15)

35

Since the displacements of the body are approximated by the shape
functions. the distributed displacement field in the u direction of an

r node element can be expressed as
Iu<x.y.z>‘°’) = I}: Nitx.y.z)ui}(°) = {INIT{u}}(°) . (A.16a)
The v and w distributed displacements are transformed into discrete

displacements in the same manner u is. and the resulting discrete

displacement field becomes

n(x.y.z) W [NTtuii W

{a} ‘0 = v(x-.y.z) [NTtvn (A.16b)

'(XIysZ) [NT{'}]

Using equations (A.14) and (A.16b) equation (A.13) becomes

{an‘°’/as} = {o} = jjj;(nir<°?ic1€°>(si<°>(.;$°?.v<e>
' IIJVIBJT(°'[CJuj‘hol‘Q'dV'e'

~jfl,tnmf1‘°?av$°? -Ij,[n]{rf}f°?d81°? (A.17a)
where
{an‘°)/as}T = [8n(°)/3u an‘°)/av au‘°’/aw} . (A.17b)

Equation (A.17a) is sometimes simplified to the following form

3‘?
[k](°)(si‘°> = (F,}(°) + {FB}(°) + (FT.(o) (A.17c)
where
[k](°? - stiffness matrix at element e

{5}(°) = displacement vector at element e

{F.}(°) = initial force vector at element e

{FB}(°) = body force vector at element e
{FT}(°) 8 surface loading force vector at element e.
/\

Y.V

 

 

 

X.U

Figure A.l Two dimensional triangular element

Equation (A.17a) is best illustrated by a two dimensional example

shown in Figure A.1. The problem consists of a single triangle

37

subjected to three leading cases. To formulate the stiffness matrix.
the shape functions must first be defined. This is accomplished by
using a linear variation of displacements u and v across the element.

resulting in two equations.

u(x.y) = N1111 + Na“: + N,u, (A.18a)

v(x.y) a N1“ + mg, + N3v, (A.18b)

In addition to equations (A.18a) and (A.18b). a third condition is

imposed. requiring the shape functions to sum to one.

N1 + N: + N, a 1 (A.18c)

Solving equations (A.18a). (A.18b). and (A.18c) for N1, N,, and N,

gives the shape functions in terms of the Cartesian coordinates. They

are
N,(x.y) 8 (1/2A)(a1 + bix + c1y) (A.19a)
N,(x.y) 8 (1/2A)(az + b,x + c,y) (A.19b)
N,(x.y) 8 (l/2A)(a, + b,x + c,y) (A.19c)
where
5 :1 Ya
2A 8 Det 1 x, y, 8 2(Area of triangle 133) (A.19d)

1 xs Ys

38

and

aa 3 x1?: ' xsis: b1 3 Ya ' Va: 01 3 x3 ' xs (A.19e)
‘s ' xs71 ' XaYs: be 3 Ya ' Yin cs 3 xi ' xs (A.19f)
‘s ' 3173 ' ‘aYin be 3 Ya ’ Ys- cs ‘ xs ' x1 . (A.19g)

Applying the (A.19) equations to the triangle in Figure A.1 results in

the following shape functions:

N; 8 (1/46)( 75 - 7x - 5y) (A.20a)
N, 8 (1/46)( -1 + 5x - 3y) (A.20b)
N, 8 (1/46) (-28 + 2x + 8y) H.200)

The matrix relating the strains to the displacements from equation

(ANS) is formed using

3N./ax o aN./ax o aN./ax o
[B] 8 0 3N1/ay 0 ‘ 3N3/8y 0 8N,/3y (A.21a)

3N1/ay aN./ax aN,/ay aN,/ax aN./ay aN,/ax
which is
1 o b. o b, o

[s] = (1/2A) 0 c1 0 c, o o. (A.21b)

c1 b1 c, be c, be

39

Applying equation (A.21b) to the element shown in Figure A.l produces

the following [B] matrix

-7 0 5 0 2 0
[B] = (1/46) 0 -$ 0 -3 0 8 (A.210)

-5 -7 -3 5 8 2 .

If the element is considered in plane strain. then the material matrix

[C] is

l-u p 0
[C] 8 E/[(1+u)(1-2u)] u 1*“ 0 (A.22a)

0 0 (l-2u)/2

where E is Young's modules and u is Poisson's ratio. If E#1000 and

“8.25 then the material matrix in equation (A.22a) becomes

1200 400 0
[C] = 400 1200 0 (A.22b)
0 0 400 .

The stiffness matrix from equation (A.17a) is

In“) - “LINN”[C](°)[B](°)dv(°) (A.23a)

but since neither the [B] matrix or the [C] matrix in this case are

40

functions of the volume. the element stiffness matrix can be rewritten

[k](°) = [BJT‘°’ICJ‘°’[B]‘°?IIIEV
g [B]T(e)[c](e)[3](e)f(e)A(e) (A.23b)

where A(°) is the element area and t(°) is the element thickness,
which usually assumes a value of one for plane strain problems. Using
the matrices in equations (A.21c) and (A22.b) in equation (A.23b) and
knowing the area of the element is 23 units. the element stiffness

matrix for the element shown in Figure A.1 becomes

' 747.8 304.4 -391.3 -17.4 -356.5 -287.0'"

304.4 539.1 -17.4 43.5 -287.0 -582.6

[k](e) . ,-391.3 -17.4 365.2 -130.4 26.1 147.8 (A.23c)
V i -17.4 43.5 -130.4 226.1 147.8 -269.6

-356.5 -287.0 26.1 147.8 330.4 139.1

 

 

.r287.0 “582.6 147.8 ~269.6 139.1 852.2.“ .

The global stiffness matrix of a structure is the sum of the element

stiffness matrices. and is found by
[K] - E“ [k](°) (A.24)

where [k](°) is summed with respect to the node numbers. The initial

force vector at element (e) is (from equation (A.17a)

41
{Fo}(e) 3 II v[B]T(°)[C](°)[8,}(°)dv(°) . (A.258)
For this element. [B] and [C] are not functions of the volume. and

since the column vector of initial strains is rarely a function of the

volume. equation (A.25a) can be expressed as

{F.1‘°’ = [81T‘°’[c1‘°>{a.1‘°’j]]6v
= [81T“?[c1‘°’{6.}‘°)£(°)A‘°’ . (A.25b)

Using an initial element strain of

[8,}T = [sxo 8y. nyo}

= [.005 .004 .001} (A.25c)

and the values of [B] and [C] as defined in equations (A.21b) and

(A.22b). the initial force vector at element (e) is
{F.}T$°) = [-28.6 —19.8 17.8 -8.2 10.8 28.0] . (A.25d)
The second type of loading that the element can undergo is body force

loading. This type of loading is usually gravity. he body force

term from equation (A.17a) is
{FB}‘°) - III;INJ{33}‘°’dv<e) (A.26a)

where

 

-n—o q-h*-‘*‘-m .. .-

42

"N1 0
0 N1
' o N,

N,o

 

 

and

{F‘}T= {x‘ 1"} . (A.26c)
If the element shown in Figure A.l has a density of p = .0075 and
gravity (g8 -386.4) is acting in the y direction. then the body force

is

um'r . {o ”mum

Equation (A.26a) becomes
(1'31“) =- [Lia lug-66.7) o N,(-66.7) o N,(-—66.7) 0}Tt(°’dA‘°’ (A.260)

and transforming the shape functions back into natural coordinates the

integration for a typical entry of (A.26e) is

”‘NiCJWdAM = t(°)A(')C./3 . (A.26d)

43

Applying the equation (A.26d) to all the entries of equation (A.26e)

results in the following body force vector
{FB}T(°’ = { o -22.2 o -22.2 o -22.2} . (A.26e)

The third force vector to be evaluated is the surface loading force
vector which is caused by boundary traction. The boundary traction

vector from equation (A.17a) is
{FT](°) = I];[NJ{I‘1‘°>aS{<e)1 . (1.27.)

The most common boundary traction term is a point load. If point 2 on

the element is loaded with point loads of

{1‘} = {-20, 15.1 . (A.27b)

Then the matrix [N] becomes

 

-0 0.
o o

[N] - 1 0 (A.27c)
o 1
o o
o o_,

 

since the shape functions. by definition, have the following property

 

44

1 at node i

N- = (A.27d)

0 at node j

the surface leading force vector will be the point loads. i.e,
{FT}T(°>= {o 0 -2o. 15. o 0} . (A.27e)

If the boundary traction term is a pressure. then it is represented by
equivalent point loads on the nodes which share the face the pressure
is acting on. Once all of the force terms have been calculated for

the element. they are added together.
m“) = um“? + {8319? + {Ir-r1“? (A.28a)

A global force vector is constructed in the same manner as the global

stiffness matrix is. thus
F a (a) . .2
U. 2.1??? . M .1.)

In the global force vector. as in the global stiffness matrix. the
element vectors are summed with respect to the nodes. The final
equation relating the global stiffness matrix and the global force

vector is

[K]{5} = [F] (A.29)

45

where {5} is the column vector of all displacements of the system.

APPENDIX,B

SHAPE FUNCTIONS OF SOLID ISOPARAMETRIC ELEMENTS

The shape functions for all solid isoparametric elements are
derived using the same method. For any n node isoparametric solid

element. the displacement in any axis direction is approximated by

. z I =- 2‘ mm - mm) (8.1)
- ‘ 1:1 . . .
where a 3 actual displacement value
3 a approximated displacement value
Ni 8 shape function of the ith node

‘i = displacement of the ith node in the direction of interest.

Th0 functions N1. Nj..... and Ni are chosen to give the appropriate
nodal displacements when the coordinates of the corresponding nodes
are inserted into equation (8.1). For a three dimensional element.

the shape function will have the following property:
1 for i equal to j

Ni]’j"j'tj)_ 3 (3.2')

0 for i not equal to‘j.

46

 

47

Furthermore. the shape functions will always sum to unity at any

location in the element.

Nifr.s.t) + Nj(r.s.t) + .... + Nm(r.8.t) = 1 . (B.2b)

One way to derive the shape function. is by assuming the variation of

the shape function is dependent on r. s, and t.

3 c d e
Nj(r.s.t) 2:;1bir s t . (B.2c)

The application of equations (3.2) is best illustrated by the element

shown in Figure 8.1.

/N

Y.V

 

T NODE K
(S=1.T=1)

 
  
  
 
   

Ln

 

NUDE [
(5"15T"[)

NODE J
(S'loTa’l)

 

x.u 7
Figure 8.1 Four node two dimensional solid element

48

This element consists of a four node two dimensional solid orientated
in it's natural coordinates (all of the examples in this appendix are
located in natural coordinates. the transformation back into spatial
coordinates will be discussed later). Equation (B.2c) applied to
Figure 3.1 will produce four equations. Using the smallest values of

d and e possible. the equations are

N1 = 1:1 + b’s + b,t + b‘st (B.3a)
NJ 3 b1 + bzs + b,t + b‘st (B.3b)
Nx - b1 + bzs + b,t + b‘st (B.3c)
NL = b1 + b,s + b,t + b‘st (B.3d)

and since the nodal values in natural coordinates are either 1 or -1.
the values 0f b3 thru b. can be solved for using equation (B.2a).
These values will not necessarily be the same for the different (B.3)
equations. Solving for the b's and rearranging the terms. the shape

functions become

"I - (l/4)(1-s)(1-t) (B.3e)
NJ - (1/4)(1+s)(1-t) (B.3f)
"I = (1/4)[1+s)[1+t) (B.3g)
NL - ulna-3711”) . (3.81.)

A typical shape function for this element is shown graphically in
Figure 3.2. The next shape function to be considered is an eight node

two dimensional isoparametric solid as shown in Figure 3.3. The shape

 

49

functions are formulated in the same way as before. only this element

has eight of them.

 

 

 

 

 

Figure B.2 Shape function of a four node element

 

 

 

X.U
Figure 3.3 Eight node two dimensional solid element

50

The shape functions for the eight node solid are

2
H
II

(1/4)(1-s)(l-t)(-s-t-1)
N3 = (174)[1+s)(1-t)[ s-t-l)
a. - unsung... as.)
NL = (l/4)[1—s)[1+t)[-s+t-1)
N, - (172571-357.-.)
NN - (1]2)(1+s)[1-tz)
"o = <1/2>(1-.2$<1+t5

NP = (1/2)(l-s)(1-t2) .

(B.4a)
(B.4b)
(B.4c)
(B.4d)
[B.4e)
(B.4f)
(B.4g)

(B.4h)

A typical shape function for this element is shown graphically in

Figure B.4.

 

Figure B.4 Shape function for an eight node element

51

The three dimensional theory for the shape functions is a continuation
of the two dimensional theory. Another variable is added to the shape
function formulation. An eight node. three dimensional solid is shown

in Figure 3.5.

 

 

 

X.U

Figure 8.5 Eight node three dimensional solid element

The shape functions for the three dimensional eight node solid element

are
N; . (ll8)(l-s)(l-t)(1-r) (B.5a)
"I = (1/8)(1+s)(1-t)(l-r) (B.5b)
N‘ = (l/8)(1+s)(1+t)(l-r) (B.5c)
NL ' (1/8)(1-s)(1+t)(l—r) (B.5d)

Nu - (1/8) (1-s)(l-t)(l+r) (B.5e)

52

NN = (1/8)(1+s)(1-t)(1+r) (8.58)
"o = (1I8)(1+s)(1+t)(1+r) (3.58)
Np = (1/8)(l-s)(l+t)(1+r) . (3.51:)

A twenty node. three dimensional solid is shown in Figure 3.6.

 

 

 

 

 

 

X.U
Figure 3.6 Twenty node three dimensional solid element

The shape functions for the twenty node three dimensional solid

element are

"I = (1/8)(l-s)(1-t)(l-r)(-s-t-r-2) (3.6a)
NJ 8 (l/8)(l+s)(1-t)(l-r)( s-t-r-2) (B.6b)

"x - (1/8)(1+s)(1+t)(1-r)( s+t-r-2) (B.6c)

53

FL = (1/8)(1-s)(l+t)(1-r)(-s+t-r-2) (B.6d)
NM = (1/8)(l-s)(1-t)(1+r)(-s-t+r-2) [B.5e)
N" = (1/8)(1+s)(1-t)(1+r)( s-t+r-2) (3.6:)
N0 = (1/8)(1+s)[1+t)(1+r)( s+t+r-2) (3.6g)
”p = (1/8)(1-s)(1+t)(1+r)(—s+m—2) (B.6h)
NQ = (1/4)(1-.2)(1-c)(1-:) (B.6i)
NR = (1/4)(l+s)(l-t2)(1-r) (B.6k)
N. = ..)....-.2;..+.,;.-.s
NT . (1)4)(1-s)(1-t2)(1-:) (B.6m)
Nu = (ll4)(l-s2)(1-t)(1+r) (B.6n)
NV 3 (l/4)(1+s)(1-t2)(1+r) (B.6o)
N, . (1/4)(1-.2)(1+:)(1+:) (B.6p)
Nx - (l/4)(l-s)(1-t2)(l+r) (8.68)
N, - (1/4)(1-s)(1-t)(1-r2) (B.6r)
NZ - (1/4)(l+s)(1-t)(1-r2) (8.6s)
NA - [174)[l+s)[1+t)[1-r2) (8.61:)
"8 - (1/4)(1-s)[1+t)[1-r2) . (B.6u)

Condensing the above elements to new elements with less nodes is based

on the following equation

a1. = (1/2)(acl + acz) (8.7a)

where n represents r. s. t. u. v. or w. the subscript m denotes a
midside node. and the subscripts c1 and 02 denote the nodes that are

collinear with node n. An eight node solid condensed into a five node

54

solid is shown in Figure 3.7.

 

 

 

XAJ

Figure 3.7 Five node two dimensional solid element

The displacement in the u direction for the eight node solid is

The five node element has three condensed sides. which will produce

three equations from equation (B.7a) in the u direction.

nu '3 (1,2) (“I 4' Ex) (3.70)
“o - (1/2)(ux + uL) (B.7d)

up - (1/2)(u.L + u1) (3.70)

55

Inserting equations (B.7e). (B.7d). and (B.7e) into equation (B.7b)

and rearranging the terms yields

u = [N1 + (Np/2)}:I + [NJ + (NN/2)]u1 + [NK + (NN/z) + (No/2)]ux +

[NL + (No/2) + (NP/2)}.L + um, (B.7f)

and the shape functions for the five node element becomes

NI --- N1 + (NP/2) (3-78)
N5 .-. N, + (Nu/2) (B.7h)
"1': = NI + (Nu/2) + (No/2) (8.71)
“f. = NL + (No/2) + (Np/2) (B.7j)
N," .. Nu . ' ‘ (B.7k)

The new shape functions derived in the v direction are exactly the
same.

If the displacements within an element are known. the strains at
any point can be determined. The displacements and strains are

related by the equation
{8} = [L]{a} (3.8:)

where [L] is a linear operator. Since the displacements are

approximated by equation (B.1). equation (B.8a) becomes

[a] 2 [lelm (B.8b)

56

where

[B] = [L][N] .

(B.8c)

The entries of matrix [B] are called the strain shape functions. The

linear operator [L] varies depending on if the problem is two or three

dimensional. In a two dimensional problem. [L] is defined as

alas O
[L] ' 0 a/at .
a/at alas

At node i. equation (Bs8c) becomes
[Bi] 3 [Lqu
or

381/38 0
[Bi] 3 0 aNi/at e

aNi/at aNiZOs
For a p node element. the [B] matrix is

[81 = [[811 raj] .... [Bpl] .

(B.8d)

(B.8e)

(3.8 f)

(B.8g)

57

If the finite element formulation is a three dimensional problem. the

linear operator [L] is

“alar 0 0 '
0 a/as 0
[L] = 0 0 an): . (3.811)

a/a: ale. 0

0 a/as a/at

 

 

_a/a: 0 aldt-

At node i. the strain shape function matrix is

 

 

’aNi/a: 0 0 ‘
0 aNi/as 0
[81] - o o aNi/at . (B.8 1)
aNi/a. aNi/a: o
o aNi/dt aNi/as
.aNi/at 0 aNi/a:.

For a three dimensional problem. the element [B] matrix assembles
according to equation (B.8g).

The above elements have all been defined in natural coordinates.
They are frequently called parent elements. Transformation of all
isoparametric elements from spatial domain to natural domain occur
using the shape functions. i.e.. the shape functions defining

displacements will also be the same function used to map the element

58

from a spatial coordinate system to a natural coordinate system. The

element shown in Figure B.8 is mapped into the parent element.

 

 

 

 

 

 

 

 

 

Figure 8.8 Four node element mapped into natural coordinates

This usually consists of a change of area. and a rotation. The

equations defining this mapping are

I = a1 + .3; + ht + a.“ (B.9a)

Y = 171 + 1),: + b,t + b‘st (3.91))
with the boundary conditions of

X(S.t) " X(-1.-1) = xI a: 4

'X( 1.

X(-1.

Y(s.t)

Y( 1.

Y("15

Solving equations

£1 8 (1/4)( :1
C, 8 [1I4)(-x1
u=iuan
a‘ 8 (1/4)( x1
b1 8 (1/4)( yI
b3 - (1/4)(-y1
b3 - (1/4)(-y1

b4 8 (1/4)[ yI

0!

X( 1,-1) 8-

1)8

1)8

Y(-1.-l) 8

Y( 1,-1) 8

l)8

1) =

H
a.
ll

so

59

(3.9 a) and (8.91))

- XI

X(s.t) 8 (ll4)(12 -

Y(s.t) 8 (1/4)(12 -

+xx+
+xx-
+xx+
.q-
+yx+
+yx-
+yx+

xL)
xi)
:1)
xx)
FL)
7L)
571,)

71..)

6s - 2st)

6t 4- 2st)

results in

12/4
«A
6
-2/4

12/4

-6/4

2/4

(B.9c)
(B.9d)
[B.9e)
(B.9 f)
(B.9g)
(8.9 h)
(B.9 i)

(B.9J)

(B.9k)

(3.91)

Once the mapping has been defined. the relationship between the strain

shape function in natural coordinates and the strain shape function in

60

spatial coordinates can be established. Using the rules of

differentiation. the strain shape functions at node B become

ENS/as = (aNfl/ax)(3x/as) + (3NB/ay)(3y/as)

aNB/at 8 (BNB/ax)(3x/3t) + (BNB/ay)(3y/3t)

or in matrix form

{Mia/6:} [ax/3s ay/as] {aura/ax}
ONB/a ax/at Bylat aNB/ay -
The above matrix equation is also written as

BNB/as aNBIGx

. = [:1

aNfl/at BNB/ay

partial

(B.10a)

(B.10b)

(B.10c)

(B.10d)

where [J] is the Jacobian matrix. To find the spatial derivatives.

the Jacobian matrix is inverted.

ORB/ax aNB/as

- [,fl .

ass/a, 3813/81;

The mapping in three dimensional problems (is similar

dimensional problems. equation (B.10c) becomes

(B.10e)

to two

 

61

3NB/as ax/as Bylas az/as ONB/ax
ONB/at = 3x/at Bylat az/at BNB/ay (B.10f)
aNB/ar axlar ay/Br az/ar aNB/Bz

or
BNB/as aNfl/ax
aNfl/at; = l J ] BNB/ay‘ . (3.108)
BNB/ar aNfl/az

The Jacobian matrix is inverted as in equation (B.10e). and the result

is
ans/ax aNB(as
'1
aNB7ay = [ J ] ans/at (B.10h)
aNB/az GNB/ar

where the Jacobian matrix (and it’s inverse) are 313 in size. With
the transformation from spatial coordinates to natural coordinates
defined. the element stiffness matrices can be transformed into
natural coordinates. the benefit being easier numerical integration of

the element stiffness matrix.

APPENDIX C

NUMERICAL INTEGRATION USING THE GAUSS-

LEGENURE.QUADRATURE METHOD

The stiffness matrix in most finite element routines is
numerically integrated by a technique called the Gauss-Legendre
Quadrature method. This method estimates the value of the integral.
f(x) by approximating the actual function with an nth-degree

interpolating polynomial pn(x) and integrating.

- d + lb 6 (C.1)
£f(x)dx Epnk) x .En(x) x

The error term is En(x) for the numerical integration. The

interpolating polynomial is of Lagrangian form. i.e.

pn(x) = f[x.] + (xrx.)f[x1.x.] + (x-x.)(x-x1)f[x,.x1.x.]

+ °"°°" + (x-x0)(x-x1)eeeee(x-xn.1)f[xnlo001:1DXo] (Cez‘)
or. expanding the first three terms

Pn(x) = f(x.) + [(:-:.)/(:.-:.)]t(:.) + [(x-x.)/(x,-x.)]r(x,) +

[(x-x.)(x-x,)/(x.-x.)(x.-x,)]£(x.) + [(x-x.)(x-x.)/(x,-x.)(x.-x.)]f<x.)

62

63

+ [(x-x.)(x-x1)/(x,-x.)(x,-x,)]f(x,) + .. ......

(C.2b)

The terms in equation (C.2b) can be collected and the equation can be

rewritten as

Pn(x) a [ Tngo'j*o[(x-xj)l(x.-xj)]]f(x.)
n
+ [”j:o'j*1[l(lx-xj)l(x1~xj)]]f(x,_)
+ + [HF0 j“I(:;-:.:J.)/(::,,-:.J.)flux“)

Equation (C.2c) is sometimes compressed to

Pn(x) 8 §:‘o[‘rT:;o’j#1[[x-xj)((xi-xj)]]f(xi)

or
n
Pn[x) - 21'0L1[x) f[xi)
where

L n
in) . n1=0.1*i[(x-xj)/(xi-xj)]

The error term from equation (C.1) is

n
En(x) 8 [.Tri‘ofx-xi)]f[x.xn.xn_1....

or

(C.2c)

(C.2e)

(C.2f)

(C.2g)

(C.3a)

64
n

Enu) = [TTigo(x-xi)]f(n+1)(n)((mtl)! (C.3b)
where u must lie on or within the limits of integration. The value of
n. however. is unknown. The next step is to transform the function
and the limits of integration into a natural coordinate system. This
is accomplished by defining a new coordinate system r as

r 8 [2x-(e+b)]/(b-a) . (C.4a)
The new function in natural coordinates will be

F(r) 8 f(x) 8 f((l/2)((b-a)r+(a+b)) . (C.4b)

Equation (C.1) in natural coordinates becomes

F(r) . 2: 01'1““ +[T]': 80(r-ri)]Fn+1[fi)[(n+l)! (C.5a)
where

Li(r) 8 W310 jii [(r- --r j)/(ri-rj )] (0.51))
and

-1 < 'i <1 . (C.5c)

T30 ti's are point transformations of the 11's. Assuming f(x) is a

65
polynomial of degree 2n&1. then 23=0Li(r)F(ri) is a polynomial of
degree n at most. Furthermore.‘TT2go(r-ri) is a polynomial of degree
n+1. thus forcing the term F(n+1)(fi)/(n+1)! to a polynomial of degree
n. This is represented by
F(“l’(fi)/(n+1)z = gn(r) (c.6)
where gn(r) is a polynomial of degree n. Equation (C.5a) becomes

F(r) 8 2:;0Li(r)F(ri) + [IrT:;o(r-ri)]gn(r) (C.7)

Integration of both sides of (C.7) yields

2““ ' ]:1§:.0L1‘r>F<=i>dr + I:1[n:go(r'ri)]xn(r)dr . (C.8)

Dropping the error term and taking the summation operator outside of

the integral produces

11%)"? F( )jlum
-1 r " 1-0 ti -1ir r

1: 280'1F(ri) (C698)
where
w [1 1. ( )dr [1 TTIII [( )/( - )]d (c 9b)
1 _1 i r -1 iso’jfii r rj V Ii Ij ‘ r. e

The term 'i is frequently called the weighing function. The location

66

of the numerical integration points is determined by the error

function. The error function is

Ii1[‘TT:;o[r-ri)]gn(r)dr . (C.10)

The best way to make the error function vanish is to select values of
‘1 that drive the function to zero. The first step is to expand the
two Polynomials (gn(r) and 'TT?.1(r-ri)) into Legendre polynomials.

The first part of the error term becomes

n
11-1'0(r-ri) 8 c.Q.(r) + c101(r) +

‘ n+1
.... + chn(r) + cn+‘03+:(r) 8 EigociQi(r) (C.11a)

and the second part of the error term becomes

3n(r) 8 d.Q.(r) + d101(r) + ... + ann(r)

' . (C.11b
2-041%“) )

Substituting equations (C.11a) and (C.11b) into the error term.

equation (C.10) results in

fi1[§:.oz.o°idjqi(r)qj(r) + cmizsodififirmmdflkx . ((3.12)

Since the Legendre polynomials have the following orthogonal

preperties

67
1
I;10n(r)qm(r)dr = 0 n 8 m (C.13a)

1
Llonumnum as o n = m (C.13b)

all terms of equation (C.12) that are of the form

1
cidleoiumJ-(nar 1 7‘ j (C.13c)

will be zero. The error term will then become

Ii1[‘rr:;o(r-ri)]gn(r)dr 8 §:=ocidiI:1[Qi(r)]2dr (C.14)

One way to make equation (C.14) equal zero is to assign a zero value
to the first n+1 01's. The coefficient cn+1 is not known. but

equation (C.11a) produces

11
TT 0(r-ri) = swank) . (c.1s)

is

In order for the left hand side of equation (C.15) to be zero. r1 must
be the roots of the Legendre polynomial. The general recursion

formula for Legendre polynomials is

Qn(r) 8 [(2n-1)/n]r0n_1(r) - [(n-l)/n]0n_,(r) (C.16a)

with

 

68

Gl.(r) 8 1 (C.16b)

Q,(r) 8 r . (C.16c)

lost finite element routines use a three point integration scheme. and

the Legendre polynomial corresponding to n83 is

93(2) = (1/2)(5r3 — 3:) . (C.16d)
The roots of equation (C.16d) are

r, = -0.77459 66692 41483 .
r, - 0.00000 00000 00000

r, 8 0.77459 66692 41483
and the weight factors from equation (C.9b) are

'1 8 0.55555 55555 55556
'3 8 0.88888 88888 88889
w

3 8 0.55555 55555 55556

The transformation from one dimensional functions to three dimensional
functions is an extension of one dimensional theory. The function is
evaluated in the other two directions in the same manner it was

evaluated in the first direction.

[i1ji1Ii1f]r"'t?drd’dt 7'2:.1§:.1§:.1'i;1;k3(r1"j(‘1’ (C-17)

69

In equation (C.17) the same number of integration points don't have to
be used in each direction. but for most three dimensional elements. an
integration grid of 31313 is used. When the stiffness matrix is
transformed from spatial to natural coordinates. a constant term is
multiplied into the entire stiffness matrix to preserve equality. The

term

IZIYI;[B(x.y.z)]T[C][B(x.y.z)]dxdydz (0.18)

will be transformed into natural coordinates as

1 1
I11] 1! 1[B(r.s.t)]T[C][B(r.s.t)]Det[J]drdsdt (C.19)

where DetIJ] is the determinant of the Jacobian matrix.

APPENDIX D

WAVE FRONT'SOLUTTON

The ANSYS finite element program uses a wave front solution
procedure.' This is a solution technique for a system of simultaneous
linear equations derived by the finite element method. The wave front
at any point is the number of equations which are active at that
instant. The ordering of elements is critical to minimize the' wave
front. for reasons of efficiency and problem size. The node numbering

is arbitrary. The active equations are

EE‘IKkjuj 3 Pk (D.1)
where xkj 8 stiffness term kj

uj 8 nodal displacement j

Fk 8 nodal force k

k 8 row number

j 8 column number

L 8 number of equations.

The elimination of an equation (i8k) begins by normalizing the first

equation in the following manner

70

71
L
§j=1(xij/Kii)uj = (Pi/xii) . (0.2)

If the finite element formulation is correct. the diagonal entry xii

will never be zero. Equation (D.2) is rewritten as

L 0 o -

Ejslxijuj = pi (0.3)
where

Kij ' Kij/Kii (B.4)
and

F1 8 171/311 . (0.5)

Equation (0.3) is saved for a back substitution solution. The

remaining equations are modified such that

a o
Kkj 3 xkj - In!“ (D . 6)

where k 8 i. These equations are assembled as '

L o o
8 F D.8
§j_zxkjuj 1 ( )

where k.varies from 2 to L. Once the first row is eliminated. the

72

other rows are eliminated by repeating the process. An example of
this solution technique is shown below. using a three by three system

of equations. The equations are

kii kis kis “1 f1

ksi tax has “a ‘ f3 . (D.9)
ksi ks: kss “a f:

The first equation is normalized so that

“i + (his/kin)“: + (k13/k11)u3 = fi/kii (D.103)

The normalized.components are then

k;, a k13/k11 (D.10b)
kis ' kll/k11 (D.10c)
f; - £,/k11 . (D.10d)

The second iteration using the wave front solution produces the

following components:

ké. = [1,, - k..<k../k..)1 (D.11a)
his - [has - k31<k1./k..)1 (D.11b)
k]. - [k,, - k,,(k,./k,.)] (D.11c)
k1. - [1,. - 1,1(k,./k,.)1 (D.11d>

f; ' [fa ‘ ksi(fs/kii)] (D.110)

73
f; = [f5 ’ k31(f3/k11)]

which are assembled as

[14. kéj {as} {:2}
his k; ns t: -
The first equation of (D.11g) is modified to become

na + [[kss ' k21(kss/k11H/[kss ’ k31(k11/k11)]]n3

a [£3 ’ k31<f1lk11)]/[k33 - k31(k13/k11)]

while the second equation becomes

stns ' f:

The term u, is solved for as
us 3 fg/k:
where

f: I [[k33 ’ k::(k13/k1i)][fg - k,1(f;(k11)]]

' [[kss ' k81(k13/k11)][f3 ' k21(£1/k11)]]

and

(D.11f)

(D.113)

(D.12a)

(D.12b)

(D.12c)

(D.12d)

 

74

k2. = [[kss ’ k31(k13’k11)][k81 ’ k31(k13/k11)]]

- [[kgg ’ kgg(kgglkgg)][kgg ’ k31(kgg/k11)]] (”.120)

The result is the same that would be obtained using some other method
such as Cramer's Rule. The advantage of the wave front technique is
the optimization of computer time during the equation solution phase

of the finite element analysis.

 

   

APPENDIX E

PLANE STRAIN

Plane strain is a specialization of three dimensional linear
elastic theory. It represents the situation where the component of
displacement normal to the x-y plane is zero. The material stiffness

matrix for the three dimensional case is

 

 

~1-p p u 0 0 0 ‘
u 1-u p 0 0 0
[C] = E/{(1+n)(1-Zu)} n u 1-n 0 0 0 - (E-l)
. [.40 ". I. 0 0 0 (1-28)/2 0 0
0 0 0 i 0 (1-2u)/2 0
L.o 0 0 0 L 0 2(1-2p)/2.

Using the fact that there are no out of plane displacements. the

material stiffness matrix compresses to

_88 u 0
[C] 8 E/[(1+u)(1-2p)] u 1-p 0 (E.2)
0 O (1-2n)/2

75

 

76

where [C] is expressed in Hooke's law as

[a] 8 [C][s] (E.3a)
where

{a)T= [ox 6y rxy] (B.3b)
and

{slT 8 [81 By Try] . (E.3c)

 

APPENDIX F

TWO DIMENSIONAL DISPLACEMENT DATA

The following tables correspond to the two dimensional model of

the wedge area from Chapter 3. They are the displacements for the

loaded nodes of the Brooks Supervillinova

modifications of various size holes.

Table F.l Y direction displacements of

Node 380.05
1 8.2990

2 8.2909

3 8.2930

4 8.2888

5 8.2870
6 8.2851'

7 8.2825

8 8.2785

9 8.2692
10 8.2427
Table F.2 I
Node 380.05
1 8.2997

2 8.2929

3 8.2959

4 8.2909

5 8.2877

380.10

8.2984
8.2911
8.2952
8.2954
8.2977
8.2953
8.2882
8.2804
8.2696
8.2427

RPO.10

-.3013
8.3075
8.3043
’.3003

380.15

8.2973
8.2917
8.2997
8.3086
8.3180
8.3148
8.2995
8.2840
8.2704
8.2425

380.15

”.3049
-.3301
‘.3308
8.3274

77

380.20

8.2957
8.2919
8.3042
8.3269
8.3569
8.3515
8.3155
8.2876
8.2709
8.2420

direction displacements of

380.20

8.3088
8.3321
8.3736
8.3771
8.3659

380.25

8.2937
8.2933
8.3128
8.3542
8.4075
8.3994
8.3386
8.2944
8.2725
8.2417

380.25

-.3170
8.3622
8.4336
8.4445
-.4268

including the

380.30

8.2914
8.2953
8.3246
-.3916
-04798
-0‘679
8.3701
8.3036
8.2748
8.2414

R-OO3O

8.3297
8.4068
8.5253
8.5502
8.5219

geometry

model with one hole

model with two holes

78

Table F.2 (cont'd)

6 8.2858 8.2988 8.3239 8.3685 8.4328 8.5324
7 8.2840 8.2946 8.3152 8.3529 8.4050 8.4829
8 8.2797 8.2853 8.2964 8.3122 8.3360 8.3694
9 8.2696 8.2713 8.2745 8.2777 8.2842 8.2932
0 8.2427 8.2429 8.2432 8.2431 8.2440 8.2453

Table F.3 Y direction displacements of model with three holes

Node 3=0.05 3=0.10 3=0.15 380.20 380.25 3=0.30

1 8.3014 8.3084 8.3234 8.3480 8.3840 8.4383
2 8.2947 8.3066 8.3302 8.3784 8.4437 8.5450
3 8.2966 8.3103 8.3378 8.3878 8.4594 8.5726
4 8.2908 8.3043 8.3317 8.3718 8.4356 8.5365
5 8.2877 8.3012 8.3282 8.3749 8.4460 8.5605
6 8.2859 8.2993 8.3261 8.3711 8.4405 8.5526
7 8.2847 8.2978 8.3241 8.3619 8.4225 8.5183
8 8.2819 8.2947 8.3192 8.3648 8.4288 8.5281
9 8.2723 8.2825 8.3019 8.3403 8.3913 8.4680
10 8.2442 8.2490 8.2588 8.2730 8.2941 8.3237

Table F.4 X direction displacements of model with one hole

Node 380.05 380.10 3=0.15 380.20 380.25 380.30

1 .01655 .01983 .02560 .03178 .04197 .05457
2 .00590 .00972 .01669 .02463 .03737 .05351
3 .01703 .02252 .03221 .04480 .06294 .08629
4 .03164 .03718 .04820 .06638 .08989 .1210
5 .04502 .04811 .05307 .06277 .07507 .09260
6 .05578 .05409 .05202 .04549 .03880 .02901
7 .06417 .06053 .05306 .03978 .02387 .00318
8 .06960 .06639 .06091 .05330 .04330 .03089
9 .07606 .07430 .07114 .06761 .06206 .05532
10

.05793 .05687 .05521 .05373 .05114 .04833

Table F.5 X direction displacements of model with two holes

Node 380.05 380.10 380.15 380.20 380.25 380.30

H

.01776 .02535 .03837 .05591 .08309 .1211
2 .00711 .01506 .03013 .05396 .08804 .1368

 

   

 

79

Table F.5 (cont'd)

3 .01718 .02341 .03411 .05170 .07827 .1191

4 .03084 .03413 .04083 .04626 .05869 .07768
5 .04440 .04566 .04801 .05145 .05790 .06804
6 .05612 .05554 .05411 .05454 .05425 .05432
7 .06423 .06065 .05499 .04387 .02942 .00783
8 .06914 .06441 .05509 .03947 .01942 8.00773
9 .07536 .07135 .06463 .05521 .04230 .02579
10

.05765 .05564 .05217 .04836 .04206 .03441

Table F.6 X direction displacements of model with three holes

Node 380.05 380.10 380.15 380.20 380.25 380.30

1 .01760 .02484 .03794 .05867 .08965 .1348
2 .00627 .01185 .02061 .03472 .05818 .09553
3 .01617 .01953 .02573 .02892 .04087 .06010
4 .03041 .03259 .03617 .04011 .05092 .06977
5 .04452 .04641 .04947 .05532 .06650 .08599
6 .05667 .05812 .06125 .06195 .06608 .07233
7 .06549 .06617 .06755 .06961 .07253 .07652
8 .07036 .06952 .06762 .06792 .06613 .06388
9 .07551 .07193 .06694 .05549 .04140 .01922
10

.05702 .05311 .04546 .03253 .01640 8.00431

APPENDIX G

THREE DIMENSIONAL CONDENSED ELEMENT'ROUTTNE

This appendix contains the computer program to generate condensed
three dimensional element stiffness matrices. The pragram was written

and run on a PRIME 750 computer.

80

81

C

C

C-8This is the driver program. STIFlONOD. that determines the
C stiffness matrix to a condensed ten node 3D isoprametric
C element. '

C

C883EAL variables

Matrix of X values of element 8 There are 20 since
a condensed 10 node element maps 20 nodes.

1(20) Matrix of T values of element

2(20) 8 Matrix of 2 values of element

X(20)

A(20) 8 Mapping coefficient from RST * A
8(20) 8 Mapping coefficient from RST ‘ B
C(20) 8 Mapping coefficient from RST ‘ C
RST(20.20) 8 Mapping matrix (see above)

Ill!"
NHM

l033(21.21) 8 Work.matrix

MU - Poissons ratio
E‘- Youngs modulas

00000000000000!)

88INTEBER variables

ITERM - The terminal location

ITYPE 8 The terminal type

IRATE 8 The terminal baud rate

ICODE - Error flag for subroutine call

00000000000

REAL X(20).Y(20).Z(20).A(20).B(20).C(20)
REAL WORK(21.21).RST(20.20).MD
GARACI‘ERH IANS

C

C

C-8Determine the terminal type and baud rate

C

CALL TERM (ITEM. ITYPE. IRATE. ICODE)

IF (ICODE.NE.0) IRATE 8 4800

N = 20

NPl 8 21
C
C-8Input the data and material properties
C
100 CALL INPUT'(X.Y.Z,E,MU)
C. . . . . ..
C88Find the location of the condensed nodes. assuming
C they lie halfway between the endpoints in all
C three directions

82

C
CALL LOCATE (X.Y.Z)
C - .
C883uild the coefficient matrix RST where
C the mapping is dependent on RST ‘ A =
C RST ..R 8 Y . RST ‘ C 8 Z

CALL ASMBL (RST)

88Solve for A. B. and C from above.
'This uses the linear equation solver
in the Math Library on the Prime.
The calling sequence is (A.X.RST.WORK.
N.NP1.IERR). The last variable. IERR.
is an error flag that indicates if
the matrix RST is singular.

OOGOOOOOOO

CALL LINEQ (A.X.RST.WORK.N.NPI.IERR)
CALL LINEQ (B.Y.RST.WORR.N.NP1.IERR)
CALL LINEQ (C.Z.RST.WORK.N.NP1.IERR)
C
C—-Check the error messgage.
C
IF (IERR.EQ.1) WRITE (1.110)
’110 FORMAT (/.1X.'888> MATRIX IS SINGULAR (888')
c. . , .- . - .
C-86ive the user the opportunity to view the
C element. This is done by mapping the element
C from natural coordinates to spatial coordin.
C
120 WRITE (1.130)
130 FORMAT (/.1X.'Would you like to see the element? ')
. READ (1.140) IANS
140 FORMAT (A1)
IF (IANS. 3031') CALL DRELE (A. B. C. IRATE)
IF (IANS. EQ.'Y') GOTO 150 ‘
IF (IANS. EQ. 'N') GOTO 150
WRITE (1.900) IANS

GOTO 120 '
C
C8-Let the user integrate the element
C

150 WRITE (1.160)
160 FORMAT (/.1X.'Would you like to integrate this element? ')
" READ (1.140) IANS
IF (IANS. EQ.'Y') CALL INTER (A. B. C. E. MD)
IF (IANS. En.'Y') GOTO 170 ‘
IF (IANS. 30. 'N') GOTO 170
WRITE (1.900) IANS '
GOTO 150 ‘
C ,
C-8Let the user input another element if they wish.

170
180

190

900

C
C

83

WRITE (1.180)

FORMAT (/ .1x.'Would you like to input another element? ')
READ (1.140) IANS

IF (IANS. BO. 'Y' ) 6010100

IF (IANS.EQ.'N') GOTO 190

WRITE (1.900) IANS

GOTO 170

CONTINUE

FORMAT (/.IX.A1.' is not a legal response. The only'.
' legal responses are Y(Yes) and N(No).')

CELL EXIT
END

C8-This subroutine. INPUT. will read in element data

C
C
C

100
110
120

130
140

150

160

170
180

and material data concerning the ten node degraded
element.

SDBROUTINE INPUT'(X. Y. 2. E. MD)
CEARACI'ER‘l IANS '
REAL X(20).Y(20).Z(20).MD

WRITE (1.110)

FORMAT (/.1X.'Enter Youngs modulas: ')
READ (1.3.ERR8120) E ‘
GOTO 130 '

WRITE‘(1.900)

GOTO 100 '

WRITE (1.140)
FORMAT (/.1X.'Enter Poissons ratio: ')
READ (1.‘.ERR8150) MU ‘

GOTO 160'

WRITE'(1.900)

GOTO 130 ' '

DO 210 J 8 1. 10. 1

IF (J.EQ. 1) IANS 8 'I'
IF (1.80. 2) IANS 8 '1'
IF (1.30. 3) IANS 8 '3'
IF (J.EQ. 4) IANS 8 'L'
IF (J.EO. 5) IANS 8 'M'
IF (1.30. 6) IANS 8 'N'
IF (1.33. 7) IANS 8 '0'
IF (1.80. 8) IANS 8 'P'
IF (1.33. 9) IANS 8 '0'
IF (1.30.10) IANS 8 'R'

WRITE (1.180) IANS. IANS. IANS
FORMAT (I. 1X. 'Enter X('. M.’). Y(' .M .'). and Z('.A1.'): ')

190

200
210

220
230

240
250
260
270

280

290
300

310

320
330
340

900
910

990

C

84

READ (1.‘ .ERR8190) X(J) .Y(J). 2(1)
GOTO 200

WRITE (1.900)

GOTO 170

CONTINUE

CONTINUE

WRITE (1.230)

FORMAT (/.1X.'The following options are available: ')
WRITE (1.240) _ '
FORMAT (/.7X.'I 8 Input new element and material data')
WRITE (1.250)

FORMAT (7X.'V 8 View current element and material data' )
WRITE (1.260)

FORMAT (7X.'E 8 Exit input mode')

WRITE (1.270)

FORMAT (I. 1X. 'Which option would you like? ')

READ (1.280) IANS

FORMAT (A1)

IF (IANS.EQ.'I') GOTO 100

IF (IANS.EQ.'V') GOTO 290

IF (IANS.EQ.'E') GOTO'990

WRITE (1.910) IANS

GOTO 220 '

WRITE (1.300) E

FORMAT (/.1X.'Youngs modulas 8 '.F12.3)

WRITE (1.310) MU

FORMAT (1X.'Poissons ratio 8 '.F8.5)

WRITE (1.320) X(1).Y(l).Z(1) '

FORMAT (/.11.'X( 1)8 '.F12.3.4X.'Y( 1)8 '.F12.3.4X.
'Z( 1)8 '.F12.3)

D0 340 J 8 2. 10. 1

WRITE (1.330) J.X(J).J.Y(J).J.Z(J)

FORMAT (IX.'X('.12.')8 '.F12.3.4X.'Y('.I2.')8 '.F12.3.4X.
'Z('.I2.')8 '.F12.3)‘ '
CONTINUE ‘

GOTO 220

FORMAT (/.1X.'Error 8 Data incompatible with program 8 '.
'Try again')‘
FORMAT (/.1X.A1.' is not a legal response. Try again.')

CONTINUE
RETURN
END' '

C88This subroutine. LOCATE. determines the location

C
C
C

C

of the condensed nodes. Their location is important
to insure complete mapping.

SUBROUTINE LOCATE (X. Y. Z)

85

DIMENSION X(20).1(20).2(20)
C , , . . , . _ _
C88First move the 3 node to it's pr0per
C twenty node location

C
X(13) 8 X(IO)
1(13) 8 1(10)
2(13) 8 2(10)
C88Node between I and K
c .
X(IO) 8 (X(2) + X(3)) / 2.0
1(10) 8 (1(2) + 1(3)) I 2.0
2(10) 8 (2(2) + 2(3)) I 2.0
C88Node between K and L
c .
X(11) 8 (X(3) + X(4)) I 2.0
1(11) 8 (1(3) + 1(4)) / 2.0
2(11) 8 (H3) + 2(4)) / 2.0
c - -. . .1 - .. .
C88Node between L and I
C
X(12) 8 (X(4) + X(l)) I 2.0
1(12) 8 (1(4) + 1(1)) I 2.0
2(12) 8 (2(4) + 2(1)) I 2.0
c ... . . . . - . .
C88Node between N and 0
C . .
X(l4) 8 (1(6) + X(7)) I 2.0
1(14) 8 (1(6) + 1(7)) I 2.0
2(14) 8 (2(6) + 2(7)) I 2.0
C88Node between 0 and P
C .
X(IS) 8 (X(7) + X(8)) I 2.0
1(15) 8 (1(7) + 1(8)) I 2.0
2(15) 8 (2(7) + 2(8)) I 2.0
c . . . . . . .
C88Node between P and M
C .
X(16) 8 (X(8) + X(5)) I 2.0
1(16) 8 (1(8) + 1(5)) I 2.0
2(16) 8 (2(8) + 2(5)) I 2.0
C .- .. - . . - . . .
C88Node between I and M
C ,
X(17) 8 (X(1) + X(5)) I 2.0
1(17) 8 (1(1) + 1(5)) / 2.0
2(17) 8 (2(1) + 2(5)) I 2.0
C , . ,

C88Node between I and N
C

 

 

 

86

X(18) 8 (X(2) + 1(6)) I 2.0
1(18) 8 (1(2) + 1(6)) I 2.0
2(18) 8 (2(2) + 2(6)) I 2.0
c . . . . . . .
C88Node between K and 0
C
X(l9) 8 (X(3) + X(7)) I 2.0
1(19) 8 (1(3) 8 1(7)) I 2.0
2(19) 8 (2(3) + 2(7)) / 2.0
C88Node between L and P
C
X(20) 8 (X(4) + X(8)) I 2.0
1(20) 8 (1(4) + 1(8)) I 2.0
2(20) 8 (2(4) + 2(8)) / 2.0
C
RETURN
END' °
C
C88This subroutine. ASMBL. assembles the coefficient
C matrix for the mapping of the function between
C natural coordinates and spatial coordinates.
C It is built on a twenty point map and passed back to
C the driver program.
C
SUBROUTINE ASIBL (RST)
REAL RST(20.20).3(20).S(20).T(20)

C88Matrices R. S. and.T are the nodal positions
C in natural coordinates.

C
DATA (3(1).131.2O)[81.0.1.0.1.0.81.0.81.0.1.0.1.0.81.0.
1 0.0.1.0.O.O.81.0.O.0.1.0.0.O.81.0.81.0.1.0.1.0.81.0/
DATA.(T(I).I‘l.20)/81.0.81.0.1.0.1.0.81.0.81.0.1.0.1.0.
1 81.0.0.0.1.O.O.0.81.0.0.0.1.0.0.0.81.0.81.O.1.0.1.0I
DATA (X(1)91.1920)[-100s-leoa-leop-1a0.1e0)1.0,1e0.1.0,
1 81.0.81.0.81.0.81.0.1.O.1.0.1.0.1.O.0.0.0.0.0.0.O.O/
C . .
m100181. 20. 1
RST(J. 1) ' 1.0
RST(J. 2) 3 3(1)
RST(J. 3) 3 8(1)
RSTU. 4) 8 TU)
RST(J. 5) ' R(J)“2
RST(J. 6) 8 S(J)“2
RST(J. 7) 8 T(J)“2
RST(J. 8) ' S(J)‘T(J)
RSTIJ. 9) ' S(J)‘R(J)

RST(J.10) 3(1)‘T(J')

87

RST(1.11) 8 S(1)“2 ‘ T(1)
RST(1.12) 8 S(1)“2 8 3(1)
RST(1.13) 8 T(1)"2 ‘ 8(1)
RST(1.14) 8 T(1)’*2 * 3(1)
RST(1.15) 8 3(1)“2 8 8(1)
RST(1.16) 8 R(1)“2 ' T(1)
RST(1.17) 8 3(1) ' S(1)"‘T(1)
RST(1.18) 8 S(1)“2 8 3(1) ‘ T(1)
RST(1.19) 8 T(1)“2 8 3(1) * 8(1)
RST(1.20) 8 R(1)“2 8 8(1) ‘ T(1)
100 CONTINUE' ' "' ' ' "
C' ..
RETURN
END~ -
C _
C
c8411. subroutine. BREE. draws the degraded twenty
C node axis by using the mapping vectors A. B. and
C C and remapping the element back from natural
C coordinates to spatial coordinates.
C
SUBROUTINE URELE (A. B. C. IRATE)
C
REAL A(20) ..B(20) C(20)
CRARACTER‘I IANS
C _ .

C88Set the initial scale value. the initial viewpoint
C and the initial reference point.
C

SCALE 8 2.

XVP 8 10.

YVP 8 10.

ZVP 8 10.

X3 8 0.

1R 8 0.

ZR 8 0.

CALL RECOVR
CALL ENITT3 (IRATE/10)
100 CALL MIDDE '
WRITE (1.110)
110 FORMAT (/.1X.'Current graphic displays set at: ')
‘ WRITE (1.120) XVP.1VP. ZVP
120 FORMAT'(5X.'Viewpoint set at (x. y.z): '.3F7.2)
‘ WRITE (1.130) XR. ER. 23 '
130 FORMAT (5X. 'Reference point set at (x. y.z): '.3F7.2)
“ WRITE (1.140) SCALE '
140 FORMAT (5X. 'Scaling factor set at: '.F7.2)

145 WRITE (1.150)

150 FORMAT (I.1X.'The following options are available: ')
WRITE (1.160)

160 FORMAT (I.5X.'D 8 Draw the element')

C

170
180
190
200
210
220

230

240

250
260

270
280

290
300

310

88

WRITE (1.170)
FORMAT‘(SX.'V
WRITE (1.180)
FORMAT‘ (5X.'R - Change the reference point')

WRITE (1.190)

FORMAT (5X.'S - Change the scaling factor')

WRITE (1.200)

FORMAT (5X.'P - View the graphic parameters')

WRITE (1.210) ‘

FORMAT (5X.'E - Exit the element draw node')

WRITE (1.220) ‘

FORMAT (I.1X.'Which option would you like? ')

READ (1.230) IANS

FORMAT (A1)

IF (IANS.m.'D') 0010 310

IF (IANS.m.'V') 0010 250

IF (IANS.EQ.'R') G010 270

IF (IANS.m.'S') GOT!) 290

IF (IANS.m.'P') GOT!) 100

IF (IANS.m.'E') com 990

WRITE (1.240) IANS

FORMAT (I.1X.A1.' is not a legal answer. Try again.')
GOT!) 145

Change the viewpoint')

WRITE (1.260)

FORMAT (l. 11.'Enter the new viewpoint (x. y.z): ')
READ (1." .RR8250) XVP. YVP. ZVP

GOT!) 145

WRITE (1.280)

FORMAT (I. 1X. 'Enter the new reference point (x. y.z): ')
READ (1'.‘ M8270) XR. YR. ZR ‘
60m 145

WRITE (1.300)

FORMAT (I.1X.'Enter the new scaling factor: ')
READ (1.8.M8290) SCALE ‘
6010 145 '

CALL REmVR

CALL DWINDO (810. 23. 10. 23. -7 .8.7 8.)
CALL CARTVP (XVP. YVP. ZVP)

CALL LOOXAT (13. YR. ZR)

CALL ZUP

CALL MAGNFY (SCALE)

C--Draw the bottom of the element
C
C

C

and then the top

W36OK81.2. 1

IF(K..m1)R8-l.
IF(K.EQ.2)3 1.

320

330

340

350
c.

360
c. .
C

C

370

380

990

89

T 8 81.

DO 320 S 8 81.. 1.. .02
CALL FCT (X. 1.2.A.B.C.3. S. T)
IF ($.30. 81. ) CALLDDVEA3 (X.1.Z)
CALL DRAWA3 (X.1.Z) '
CONTINUE ‘

S 8 1'.

DOB3OT881.. 1.. .02
CALLFCT (X.1.Z......ABCRST)
CALL DRAWAS (X. 1. 2)

CONTINUE '

T8 1.

D0 340 S 8 1.. 81.. 8.02
CALL FCT (X.1.Z.A.B.C.R.S.T)
CALL DRAWA3 (X.1.Z)
CONTINUE

S 8 81.

DOB5OT8 1..-1. 802
CALLFCT (X.1.....ZABCRST)
CALL DRAUAB (X. 1.2)
CONTINUE

CON TINDE

C—Connect the top and bottom

D038OX81.4.1

IF (1.39.1) S 8 81.
IF (X.m.1) T 8 81.
IF (3.30.2) S 8 1.
IF (X.m.2) T 8 81.
IF (LEQJ) S 8 1.
IF (3.30.3) T 8 1.
IF (1.30.4) S 8 81.

T8 1.

IF (3.80.4)

D0 370 R 8 81.. 1.. .02

CALLFCT (X.1.2..AB.C3..ST)

IF (3'. m.81. ) CALL MVEA3 (X. 1.2)
CALL D3AWA3 (X. 1.2)

CONTINUE ‘

CONTINUE
CALL ANBDDE

6010 145
CONTINUE

RETURN ‘
END

 

 

90

C
C-8This subroutine. FCT. maps points in natural
C coordinates back into spatial coordinates.

C

SUBROUTTNEIFCT (X.1.2.A.B.C.R.S.T)
C

REEL A(20).B(20).C(20)
C _

8 A(1) + A(2)8R + A(3)8S + A(4)8T + A(5)8R8R

+ A(6)8S8S + A(7)8T‘T 4' A(8)"S‘T + A(9)8S8R +
A(10)8R8T + A(11)8S8S8T-+ A(12)‘S8S83 + A(13)8S8T8T
+‘A(14)8R8T8T*+'A(15)8R8R8S'+ A(16)8R8R8T"+"‘ '
A(17)8R8S8T + A(18)8S8S8R8T + A(19)8R8S8T8T +
A(20)8R838S8T

h8kih8h8h8

Y = 3(1) + B(2)8R + B(3)8s + B(4)8T + B(5)8R8R

+ B(6)8S8S + B(7)8'r8'r + B(8)8S8'r + B(9)8S8R +
B(lO)8R8T + B(ll)8S8S8T‘+ B(12)8S8S8R + B(13)8S8T8T
+ B(14)838T8T + 8(15)8R838S + B(16)8R8R8T +
B(l7)8R8S8T-+ B(18)‘S8S8R8T + B(19)8R8S8T8T +
B(20)83838S8T ‘

~F8F8F8F8P8

2 8 C(1) + C(2)8R + C(3)8S + C(4)8T'+ C(5)8R8R

+ C(6)8S8S +‘C(7)8T8T + C(8)8S8T"+ C(9)8S8R +
C(10)838T‘+ C(11)8S8S8T'+ C(12)‘S8S8R + C(13)8S8T8T
+ C(14)838T8T'+ C(15)8R8R8S + C(16)8R8R8T'+
C(17)8R8S8T'+ C(18)8S8S8R8T + C(19)8R8S8T8T +
C(20)8R8R8S8T

h8hik8h8h8

RETURN
END
C
C
C-8This subroutine. INTEB. integrates the element
C stiffness matrix and builds the matrix.
C The integration scheme used is a 3X3X3 Gauss
C Legendre quadrature method
C
C88Variab1e identification
C
C-8REAL variables

A(20) - Mapping coefficients 8 X direction
B(20) 8 Mapping coefficients 8 1 direction
C(20) 8 Mapping coefficients 8 2 direction

3(3) 8 3 natural coordinate integration points
8(3) - 8 natural coordinate integration points
T(3) 8 T natural coordinate integration points
WR(3) 8 Weight functions. r direction
WS(3) 8 Weight functions. s direction
WT(3) 8 Weight functions. t direction

000000000000

91

DSFRST(3.10) 8 Derivative of the shape functions with resp to RST
DSFX12(3.10) 8 Derivative of the shape functions with resp to X12

BMAT(6.30) 8 The 3 matrix used to form the stiffness matrix
BTMAT(30.6) 8 The transpose of the 8 matrix

CMAT(6.6) 8 The material matrix used to form the stiffness matrix
INV1AC(3.3) 8 The inverse of the Jacobian matrix

PART(30.30) - A partial sum of the stiffness matrix
STIF(30.30) 8 The element stiffness matrix

0000000060060

SUBROUTINE mm; (A.B. C. E. MU)

0

REAL A(20).B(20).C(20).MU

REAL R(3).S(3).T(3).WR(3).WS(3).WT(3)
REAL DSFRST(3.10) .DSFX12(3.10)

REAL BMAT(6.30).BTMAT(30.6).CMAT(6.6)
REAL INV1AC(3.3) ‘
REAL STIF(30'.30) .PART(30.30).WORX(30.6)

DATA (3(1).I81.3)l8.7745966692.0.0..7745966692!

DATA (8(1).I81.3)l8.7745966692.0.0..7745966692/

DATA (T(I).181.3)l8.7745966692.0.0..7745966692/

DATA (“(1) . 181 .3) / .5555555556 . .8888888889 .‘ .5555555556/
DATA (WS(I) . I81 .3)/ .5555555556 . .8888888889 . . 5555555556!
DATA (WT(I) . I81 .3)/.5555555556. .8888888889 . .5555555556l
DATA STIF/90080.O/'

DATA PART/90080.0/

DATA EMT/18080.0/

DATA BMAT/1808O.OI'

DATA GMAT/3680.0/'

C .

C88Start the integration loop
C .

D0 330 I
DO 320 1
DO 310 K

1. 3. 1
1. 3. 1
l. 3. 1
C . .
C88Get the derivatives of the shape functions
C in thier natural coordinates

C

CALL DERSF (DSFRST.R(I) .S(X) .T(1))

C . . . . _.
C—Determine the Jacobain matrix inverse. the Jacobian
C matrix determinant. for the numerical values of

C 3. S. and T.

C

CALL JACOB (A.B.C.3(I).S(X).T(1).DETJAC.INV1A0)
C . . . .
C88Transform the derivatives of the shape functions
C into X12 spatial coordinates.

92

C This involves using the matrix multplication
C in the PRIME math library. called MMLT
C
CALL MMLT (DSFX12.INV1AC.DSFRST.3.3.10)
c .

C88Set the B matrix from DSFXIZ for the stiffness matrix
C 'Also find the transpose of the B matrix. These two
C matrices are stored in BMAT and BTMAT
C
CALL MAKER (DSFX12.BMAT.BTMAT)
C
C
C88Get the material matrix. this is called inside
C the integration loop because some matrix routines
C destroy the product matrices.

C
CALL MAKEC (CMAT.E.MU)
C A
C—Build the stiffness matrix parts
C
CALL MIILT (WORK.BTMAT.CMAT.30.6.6)
CALL MILT (PART.WORX.BMAT.30.6.30)
C

C88Multiply the weight factors and DE'DAC by the
C individual summation part of the stiffness matrix
C
CALL MSCL (PART. PART. 30. 30. (W3(I)8WS(X)8WT(1)))
CALL MSG. (PART. PART. 30. 30. DE'HAC) ‘ ‘
C
C88Sum the partial stiffness matrix to the entire term
c .
CALL MADD (STIF.STIF.PART.30.30)
C
310 CONTINUE
320 CONTINUE
330 CONTINUE
c .
C88Output the results to a file named OUTPUT
C
OPEN (14.FILE-'WTPUT')
D0 510*I 8 1. 30. 1'
DO 500 1 8 1.30. 6
WRITE (14.8) STIF(I.1). STIF(I. 1+1). STIF(I. 1+2). STIF(I. 1+3).
1 STIF(I.1+4). STIF(I. 1+5) ‘
500 CONTINUE
510 CONTINUE

C
CLOSE (14)
C . . a .
RETURN

EVID'

93

C88This subroutine. MAKEC. assembles the material

C
C

C

C

100
110

C

SUBROUTINE MAKEC (CMAT. E. MU)

REAL CMAT(6.6).MU

CMAT(1.1)
CMAT(1.2)
CMAT(1.3)

CMAT(2.1)
CMAT(2.2)
CMAT(2.3)

CMAT(3.1)
CMAT(3.2)
CMAT(3.3)

CMAT(4.4)
CMAT(5.5)
CMAT(6.6)

D0 110 I 8
D0 100 1 8

CMAT(I.J)
CONTINUE'
CONTINUE

RETURN
END

1
1

1 8 MU

MU

MU

MD

1 - MU

m .

MU

MD

1 8 MB

(1 8 (2.08MU)
(1 8 (2.08MB)
(1 8 (2.08MU)
. 6. 1

. 6. 1

( E

matrix in a location called CMAT.

) / 2
) / 2
) / 2

.0
I0
0

C-8This subroutine. DERSF. calculates the

t

(183)
(183)
(1-T)

(18R)
(1+8)
(1-T)

(183)
(183)

O

derivatives of the shape functions in

(2.088

(S + 2.

+ T + R +
08T + R +

(S + T + 2.083 +

(2.088
(2.08T
(T _ 8

(s + 2.

C

C natural coordinates.

C

C
SUBROUTINE DERSF (DSFRST.R. S. T)
REAL DSFRST(3.10)

c . .
DNDSI = .125 8 (1—T)
DNDTI 8 .125 8 (188)
DNDRI 8 .125 8 (188)

C
DNDSI 8 .125 8 (18T)
DNDTU 8 .125 8 (183)
DNDRI 8 .125 8 (1+8)

C
DNDSI = .125 8 (1+T)
DNDTI 8 .125 8 (1+8)
DNDRI 8 .125 8 (1+8)

(1+T)

(2.0.x

- T _ R _
- s + n +
+ 2.083 +

+ T — R 8
08T — R 8
8 S 8 T +

/ ((1+MU)8(18(2.08MU)))) 8 CMAT(I.1)

 

DNDSL 8
DNDTL 8
DNDRL 8

DNDSM 8
DNDTM 8
DNDRM 8

DNDSN 8
DNDTN 8
DNDRN 8

DNDSO 8
DNDTO 8
DNDRO 8

DNDSP 8
DNDTP 8
DNDRP 8

DNDSQ 8
DNDTQ 8
DNDRQ 8

DNDSR 8
DNDTR 8
DNDRR 8

DNDSS 8
DNDTS 8
DNDRS 8

DNDST'8
DNDTT
DNDRT'8

DNDSU 8
DNDTU 8
DNDRU 8

DNDSV 8
DNDTV 8
DNDRV 8

DNDSW 8
DNDTW 8
DNDRW 8
DNDSX 8
DNDTX 8
DNDRX 8

DNDS1 8

.125
.125
.125

.125
.125
.125

.125
.125
.125

.125
.125
.125

.125
.125
.125

.25
.25
.25

.25
.25

.25

.25
.25
.25

.25
.25
.25

.25
.25

.25

.25
.25
.25

.25
.25
.25

.25
.25

.25

(1+T)
(188)

(188)

(18T)
(188)
(188)

(1—1‘)
(1+8)
(1+8)

(1+T)
(1+8)
(1+8)

(1+T)
(188)
(188)

(181‘)
(183)
(18T)

(18R)
(1+8)
(1+8)

(1+T)
(183)
(1+T)

(183)
(188)
(188)

(181')
(1+3)
(18T)

(1+R)
(1+8)
(1+8)

(1+T)
(1+3)

(1+T)

(1+3)
(188)

(188)

(181‘)

94

(18R) 8 (2.088 8 T'+ R.+
(18R) 8 (2.08T 8 8 8 R 8

(1+T) 8 (s -:'r+ 2.08R

...

(1+3) 8 (2.088 + T - R +
(1+3) 8 (S + 2.08'1‘ - R +

(18T) ‘ (2.083 8 8 8 T

(1+R) 8 (2.088 - 'r + R

(1+3) 8 (2.08'1' 8 S 8 R +

(18T) 8 (S 8 T-+ 2.083

(1+3) 8 (2.088 + T + 3
(1+3) 8 (S + 2.08T + R
(1+T) 8 (S'+ T'+ 2.083

(1+3) 8 (2.088 8 T‘8 3
(1+3) 8 (2.08T 8 S'+ R
(1+T) 8 ('1‘ 8‘s + 2.0811

(183) 8 (82.0)88

(18(8882)) 8 (81.0)
(18(8882)) 8 (81.0)

(18(T882))
(183) 8 (82.0)8T

(18(T882)) 8 (81.0)

(183) 8 (82.088)
(18(8882)) '

(18(8882)) 8 (81.0)

(1-(T882)) 8 (81.0)

(18R)‘8 (82.0)8T

(18(T882)) 8 (81.0)

(1+R) 8 (82.0)88

(18(8882)) 8 (81.0)

(18(8882))

(1-(‘1‘882H
(1+R) '8 (82.0)8'1‘
(1801882))

(1+R) 8 (82.0)88
(18(8882))
(18(8882))

(18(T8e2)) 8 (81.0)

(1+R)" (82.0)8T
(18(T882))

(18(3882)) 8 (81.0)

C

DNDTY 8
DNDRY 8

.25 .
.25 8

DNDSZ 8
DND'IZ 8
DNDRZ 8

.25 8
.25 8
.25 .

.25 8
.25 8
.25 8

”(BSA 8
DNDTA 8
DNDRA 8

.25 8
.25 8
.25 8

DNDSB 8
DND'IB 8
DNDRB 8

(188) 8
(188) 8

(1-T) 8
(1+8) 8
(1+8) 8

(1+T) 8
(1+8) 8
(1+8) 8

(1+T) 8
(1-S) 8
(1—S) 8

C88Load the derivatives into

C

DSFRSTU. 1) 8
DSFRST(2. 1)
DSFRST(3, 1)
DSFRSTU. 2)
DSFRST(2. 2)
DSFRST(3. 2)
DSFRSTU. 3)
DSFRST(2. 3)
DSFRST(3. 3)
DSFRST(1. 4)
DSFRST(2. 4)
DSFRST(3. 4)
DSFRSTU. 5)
DSFRST(2. 5)
DSFRST(3. 5)

DSFRS'NI. 6)
DSFRST(2. 6)
DSFRST(3. 6)
DSFRSTU. 7)
DSFRST(2. 7)
DSFRST(3. 7)
DSFRSTU. 8)
DSFRSTQ. 8)
DSFRST(3. 8)
DSFRSTU. 9)
DSFRS'I'(2. 9)
DSFRST(3. 9)
DSFRST(1.10)
DSFRST(2.10)
DSFRST(3.10)

RETURN
mu ‘

DNDSI +
DNDTI +
DNDRI +
DNDSI 8-
DNDTI +
DNDRI +
DNDSI 8-
DNDTI 4-
DNDRI +
DNDSI. +
DNDTI. +
DNDRL +
DNDSI +
DNDTI! +
DNDRI! +

DNDSN +
DND'IN +
DNDRN +
DNDSO +
DNDTO +
DNDRO +
DNDSP 4-
DNDTP +
DNDRP +
DNDSQ
DNM‘Q
DNDRQ
DNDSI]
DNDTU
DNDRI]

95

(1-(R..2)) . (-1.0)
(1413-8 (82.0)":

(18(R882))
(18(8882)) 8 (81.0)
(1-T)‘8 (82.0)8R

(1—(R882))
(18(8882))
(1+T) 8 (82.0)8R

(1-(R882)) 8 (81.0)
(1-(R882)) '
(1+T) 8 (82.0)8R

the nut: ix DSFRST

(DNDSI/2.0) + (DNDSY/2.0)
(DNDTI/2.0) + (DNDTY/2.0)
(DNDRT/2.0) + (DNDRY/2.0)
(DNDSR/2.0) + (DNDSZ/2.0)
(DNDTR/2.0) 8 (DNDI'Z/2.0)
(DNDRR/2.0) 8 (DNDRZ/2.0)
(DNDSR/2.0) '8 (DNDSS/2.0)
(DNDTR/2.0) 8' (DNDTS/2.0)
(DNDRR/2.0) + (DNDRS/2.0)
(DNDSS/2.0) + (DNDSI/2.0)
(DNDTS/2.0) + (DNDTI/2.0)
(DNDRS/2.0) + (DNDTI/2.0)
(DNDSX/2.0) 8’ (DNDSI/2.0)
(DNDTI/2.0) + "mm/2.0)
WWII/2.0) + (DNDRY/2.0)
(DNDSV/2.0) 8’ (DNDSZ/2.0)
(DNDTV12.0) + (DNDTZ/2.0)
(DNDRV/2.0) 8' (DNDRZ/2.0)
(DNDSV/2.0) 8' (DNDSW/2.0)
(DNDTV/2.0) 4’ (DNDTI/2.0)
(DNDRV/2.0) + (DNDRI/2.0)
(DNDSW/2.0) + (DNDSI/2.0)
(DNDTI/2.0) 8’ (DNDTI/2.0)
+

“MRI! 2 .0)

(DNDRI/2.0)

++++++

++++++

(DNDSA/2.0)
(DNDTA/2.0)
(DNDRA/2.0)
(DNDSB/2.0)
(DNDTR/2.0)
(DNDRB/2.0)

(DNDSA/2.0)
(DNDTA/2.0)
(DNDRA/2.0)
(DNDSB/2.0)
(DNDTR/2.0)
(DNDRB/2.0)

C

96

C-‘l‘hia subroutine. IAmB. determines the Jacobian
matrix of the twenty node map. It also finds
the determinant of the Jacobian and it's inverse.

C
C
C

C

C

HH

SUBBOUTINE JACOB (A.B.C.R.S.T.DE'DAC.INVIAC)
REAL A(20).B(20).C(20).INVJAC(3.3).JAC(3.3).WORK(4.6)

Dxns 8 A(3) + A(6)82.08s + A18)8T‘+ A(9)8R +
A(11)82.08S8T + 'A(12)82.08S8R + A(13)8'1‘882 + A(15)8R882 +
A(17)8R8‘r'+'A(18)82.08R8S8T‘+ A(19)8R8T882 + A(20)8T8R882
DXDT‘8 M4) + A(7)82.08'1' + A(8)8s + A(10)8R +
A(11)8S882‘+ A(13)82.08S8T + A(14)82.08R8T + A(16)8R882 +
A(17)8R8S + A(18)8R8S882-+ A(19)82.08R8S8T + A(20)8S8R882
DDR8 A(2) + A(5)82.08R + A(9)8s + A(10)8'r +
A(12)8S882‘+ A(14)8T882 + A(15)82.08R8S + A(16)82.08R8T +
A(17)888T + A(18)8'r88882 + A(19)8S8'r882 4- A(20)8R8S8T

DYDS = 3(3) + B(6)82.08s + B(8)8'r + B(9)8R +
B(11)82.08S8T + B(12)82.08S8R +'B(13)8T882 + B(15)8R882 +
B(17)8R8T'+ B(18)82.08R8S8T + B(19)8R8T882 + B(20)8T8R882
DYDT'8 3(4) +'B(7)82.08T + B(8)8S + B(10)8R +
B(11)8S882'+ B(13)82.08S8T + B(14)82.08R8‘l‘ + B(16)8R882 +
B(17)8R8S + B(18)8R8S882 + B(19)82.08R8S8‘I' + B(20)8S8R882
BID! 8 3(2) + B(5)82.08R + B(9)8S 4- B(10)8T +

B(12)8S882 + B(14)8T882 + B(15)82.08R8S + B(16)82.08R8T +
B(17)8S8T + B(18)8T8S882 + B(19)8S8T882 + B(20)8R8S8T‘

M138 8 C(3) + C(6)82.08S + C(8)8T + C(9)8R +

C(11)82. 08881‘ 4' C(12)82. 08888 + C(13)8T882 + C(15)8R882 +
C(17)8R8T + C(18)82. 08R8S8‘l' + C(19)8R8T882 + C(20)8T8R882
DZDT 8 (1(4) + C(7)82. 081' + C(8)8S + C(10)8R +

C(11)8S882 + C(13)82. 08S8T + C(14)82. 0838'! + C(16)8R882 +
C(17)8R8S + C(18)8R8S882 + C(19)82.08R8S8T + C(20)8S8R882
DZDR 8 C(2) + C(5)82.08R + C(9)88 + C(10)8T +

C(12)8S882 + C(14)8T882 + C(15)82.08R8S + C(16)82.08R8T +
C(17)8S8T + C(18)8T8S882 + C(19)8S8T882 + C(20)8R8S8T

C-8Assemb1e the Jacobian matrix

C

JAC(1.1) 8 DXDS
JAC(1.2) 8-DYDS
JAC(1.3) 8 DZDS

JAC(2.1) 8 DXDT
JAC(2.2) 8 DYDT
JAC(2.3) 8 DZDT

JAC(3.1) 8 DXDR
JAC(3.2) 8 DYDR
JAC(3.3) 8 DZDR

97

C88Use the matrix inversion routine on the PRIME.

C The calling sequence is
C call minv (OUTPUTuINPUT,ISIZE,WORX.ISIZE+1,
C ISIZE-I-ISIZEJERROR FLAG)
C .
C
N 8 3
NP1 8 4
NPN 8 6
C _
CALL MINV (INVJAC.JAC.N.WORK.NP1,NPN.IERR)
C

C88Solve for the Determinant of the Jacobian
C

+

DETIAC 8 (JAC(1.1)8JAC(2.2)8JAC(3.3))
(JAC(1,2)8JAC(2.3)8JAC(3,1))
(JAC(2.1)8JAC(3.2)8IAC(1.3))
(JAC(3.1)8JAC(2.2)8IAC(1.3))
(JAC(2.1)8JAC(1.2)8JAC(3.3))
(JAC(3.2)8JAC(2.3)8JAC(1.1))

...

F8P8F8F8F8

IF (IERR.NE.0) WRITE (1.100)

100 FORMAT (I.1X. 0......) JACOBIAN ' IS SINGULAR (888')
RETURN

END‘

88This subroutine. MAKER. builds the B matrix
(called.BMAT) in the stiffness matrix.
The transpose of B (called.BTlAT) is also
assembled and passed back in INTEG.

OOOOOOG

SUBROUTINE MAKER (DSFXIZ.BMAT.BTMAT)

O

REML DSFX!Z(3.10).BMAT(6.30).BTNAT(30.6)
C88Build the B matrix from the matrix DSFXIZ
C ,

K 8 1
00100181.30.3
BMAT(1. I) 8 DSFXIZ(1.K)
BMAT(2.I+1) 8 DSFX!Z(2.K)
BMAT(3.I+2) 8 DSFXYZ(3.R)
BIAT(4. I) 8 DSFXYZ(2.K)
BMAI(4.I+1) 8 DSFXYZ(1.K)
BNAT(5.I+1) 8 DSFXTZ(3.K)
BMAT(5.I+2) 8 DSFX!Z(2.K)
BIAI(6. I) 8 DSFIIZ(3.K)
BNAT(6.I+2) 8 DSFX!Z(1.K)
I 8 K'+ l

100 CONTINUE

98

C88Find the transpose of the B matrix and place it
C in the location BTMAT
C
no 120 I 8 l. 30. 1
D011018 1. ‘6. 1
BTNAT(I.J) 8 BNAT(J.I)
110 CONTINUE ' ’
120 CONTINUE

RETURN
END ‘

 

 

 

 

“WRuin?