DYRAMIC ANALYSES OF
NQNLENEAR SPACE FRAMES

Thesis for the fiegree 9? Ph. D.
MECREQAN STATE UNE‘JERSWY
FERWEGGW FARHGOMAED
197a

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MW \\ NT

3 gas 10758 7960‘ _
WM. . Michiga‘ ltate
Ulll" 5333?] (,4

. t‘
P" - 1" r” ? u s
“I.

 

 

 

 

 

 

 

 

 

W W i, _ 12%;; ,,

‘

—I

This is to certify that the

thesis entitled

DYNAMIC ANALYSIS OF
NONLINEAR SPACE FRAMES

presented by

Fereydoon Farhoomand

has been accepted towards fulfillment
of the requirements for

Ph . D degree in Lililjngineer ing

 

Q. KW

Major professor

 

Date May 20‘ 1970

 

0-169

 
 
   
     
 

  

. . 4 .. ”I A;

r; :3? HNSIG av
“ "MG & SHIS'

3'. I "UK SINK" INC. 1;

LIL LIBRARY muons Q1

 

 

. in

x
" o‘v’

Lia‘AApuAAans m‘

. 4:?
gigw~67 .:.__.5

4!" 7‘, 010 "

W'W

I .l i' - f“-
.

, {3‘8 0‘3} 0 2

 

 

 

 

 

 

 

 

 

ABSTRACT

DYNAMIC ANALYSIS OF
NONLINEAR SPACE FRAMES

BV

.-

Fereydoon Farhoomand

In this thesis, a matrix formulation is presented
for the analysis of dynamically loaded space frames. The
effects of both material and geometric nonlinearities are
taken into consideration. These effects are restricted,
respectively, to the case of linearly elastic-perfectly
plastic materials and the case of small rotati01s, i.e.,
the case in which the rotation angles are negligible with
respect to unity.

The analysis begins with specifying two yield
condition equations for a typical member cross-section.
Then, the incremental force-displacement relations for
a space-frame member are derived for several cases.
Firstly, these relations are summarized for a linearly
elastic member with the effects of geometric nonlinearities
taken into account. Secondly, the relations for a member
whose end cross-sections are yielding are derived again

with the latter effects accounted for. Finally, the

Fereydoon Farhoomand

incremental force-distortion relations are derived for an
elastoplastic discrete model under geometrically linear
conditions.

The analysis is further carried on by describing a
mass lumping procedure that considers rotary inertia.
Then, the equation of motion for a typical "free" joint
is derived. Following that, the criterion by which a mem-
ber cross-section is ruled to be yielding is described.

It is shown next that the force vector acting on a
yielding cross-section is prevented from proceeding beyond
the yield surface. Finally, the steps of the numerical
procedure employed to determine transient response are
described. I

A computer program is prepared for the implementa-
tion of the analysis. Three numerical problems are con-
sidered: a cantilever beam, a six-member space frame, and
a two-bay two-story building frame. Dynamic loading is
either provided by concentrated loads applied to free
joints, or generated by ground motions due to earthquake.

Several comparative studies on the numerical prob-
lems mentioned above are presented. These studies show
that good agreements exist between the results provided
by the geometrically linear formulation and those given
in a published report in which a different method was

used. This may be construed as evidence for the

Fereydoon Farhoomand

validity of the present analysis. The method given here,
however, requires significantly less computer time.

The comparative studies also show that plastic
displacements as predicted by the geometrically nonlinear
formulation are generally larger than those resulting
from the geometrically linear version. But, on the whole,
when axial loads are small (as compared to the corre-
sponding Euler loads), the influence of geometric non—
linearities on numerical results does not seem significant.
However, as axial loads increase, such influence rapidly

grows.

DYNAMIC ANALYSIS OF
NONLINEAR SPACE FRAMES

BY

\
l"

Fereydoon Farhoomand

A THESIS

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

DOCTOR OF PHILOSOPHY
Department of Civil Engineering

1970

ACKNOWLEDGMENTS

The writer is deeply indebted to his academic
advisor, Dr. Robert K. Wen, Professor of Civil Engineer—
ing and Acting Chairman of the department, for his con-
tinuous advice and assistance during the past five years.

Thanks are expressed to the other members of the
writer's guidance committee, Dr. Charles E. Cutts,

Dr. William A. Bradley, and Dr. Chi Y. Lo, for their
help. The encouragement of Dr. James L. Lubkin and Dr.
Lawrence E. Malvern, who formerly served on the writer's
guidance committee is also acknowledged.

Sincere appreciation is extended to the National
Science Foundation and Division of Engineering Research

for their financial support of this investigation.

ii

TABLE OF CONTENTS

Page
ACKNOWLEDGMENTS . . . . . . . . . . . . . . . . . ii
LIST OF TABLES . . . . . . . . . . . . . . . . . . v
LIST OF FIGURES . . . . . . . . . . . . . . . . . Vi
LIST OF SYMBOLS O O O O O O O O O O I O O O O O O Vii
CHAPTER I: INTRODUCTION . . . . . . . . . . . . . 1
1.1 Objective . . . . . . . . . . . . . . . . l
1.2 Previous Studies . . . . . . . . . . . . 2
1.3 Assumptions and Limitations . . . . . . . 4
1.4 Outline of Present Study . . . . . . . . 6
1.5 General Definitions . . . . . . . . . . . 7
CHAPTER II: FORCE DEFORMATION RELATIONS
FOR A MEMBER . . . . . . . . . . . . 11
2.1 Yield Condition Equations . . . . . . . . 11
2.2 Incremental Force-Displacement
Relations for a Geometrically
Nonlinear Member . . . . . . . . . . . . 14
2.2.1 Elastic Member . . . . . . . . . . 15
2.2.2 Elastoplastic Member . . . . . . . 21
2.2.3 Elastic Return . . . . . . . . . . 26
2.3 Incremental Force—Distortion
Relations for a Discrete Model
of a Geometrically Linear Member . . . . 27
CHAPTER III: DYNAMIC ANALYSIS AND SOLUTION . . . 33
3.1 Dynamic Analysis . . . . . . . . . . . . 33
3.1.1 Joint Mass Matrix . . . . . . . . 34
3.1.2 Equation of Motion for a
Joint Rigid Body . . . . . . . . . 37
3.2 Dynamic Solution . . . . . . . . . . . . 40
3.2.1 Insertion of a Plastic
Hinge at a Cross-Section . . . . . 40
3.2.2 Force Containment at a
Plastic Hinge . . . . . . . . . . 41
3.2.3 Numerical Procedure
of Solution . . . . . . . . . . . 42
3.2.4 Choice of Time Increment . . . . . 45
3.3 Computer Program . . . . . . . . . . . . 46

iii

 

fl!'\ R_-'

— _, .
V‘s-‘8- ‘ .4

J-

hut-.J-

norms

_‘~.|u

 

 

 

Page

CHAPTER IV: APPLICATIONS . . . . . . . . . . . . 47
4.1 Cantilever Beam . . . . . . . . . . . . . 48
4.2 Triangular Frame . . . . . . . . . . . . 52
4.3 Building Frame . . . . . . . . . . . . . 53
4.4 Computation Time . . . . . . . . . . . . 55
CHAPTER V: SUMMARY AND CONCLUSIONS . . . . . . . 56
TABLE . . . . . . . . . . . . . . . . . . . . . . 59
FIGURES . . . . . . . . . . . . . . . . . . . . . 60
LIST OF REFERENCES . . . . . . . . . . . . . . . . 79
APPENDIX I: STIFFNESS MATRICES . . . . . . . . . 81

APPENDIX II: INCREMENTAL FORCE-DISPLACEMENT
RELATIONS FOR A GEOMETRICALLY LINEAR
MEMBER WITH END RELEASES . . . . . . 87

APPENDIX III: COMPUTER PROGRAM . . . . . . . . . 91

iv

LIST OF TABLES
Table Page

4-1 Comparison of Computer Time between
Reference 13 and Present Work . . . . . . 59

LIST OF FIGURES

Typical Space Frame and Notation
Different Coordinate Systems . .

End Forces and Displacements
of a Typical Member . . . . . .

Free-Body of a Differential Element . . .

Geometrical Interpretation of Flow Law . .

Discrete Model as Used in Reference 13 . .

Rigid Body Associated with Joint
Treatment of Force Vector during
Cantilever Beam and Loading . .

Transverse Displacement of
Cantilever Beam versus Time . .

Effects of Axial Load . . . . .

Axial Displacement of Cantilever
Beam versus Time . . . . . . . .

Triangular Frame and Loading . .

Response of Triangular Frame
(Elliptic Yield Function) . . .

Response of Triangular Frame
(Parabolic Yield Function) . . .

Comparison between Geometrically Linear and

J O O O O

Yielding

Nonlinear Responses of Triangular Frame .

Building Frame . . . . . . . . .

Plan View of Distorted Frame Predicted
by Geometrically Linear Formulation . .

Plan View of Distorted Frame Predicted
by Geometrically Nonlinear Formulation . .

vi

Page
60

61

62
63
64
65
66
67

68

69

7O

71

72

73

74

75

76

77

78

LIST OF SYMBOLS

Flexibility matrix;
cross-sectionsl area;
effective areas of shear;
flow-constant matrix;
flow constant;

displacement transformation matrix
defined by Equation (1-3);

modulus of elasticity;

matrix defined by Equation (2-31)
or (2-51);

axial strain;
superscript denoting "elastic part";
member force (vector);

member force (vector) due to
geometrically nonlinear terms;

axial force;

fully plastic axial force;

Shearing forces;

modulus of shear;

matrix defined by Equation (2-37);
matrices defined by Equations (2—26);
matrix defined by Equation (1—2);

third order identity matrix;

vii

Z Z

2I

branch inertia tensor expressed at
branch mass center b;

branch inertia tensor expressed at
joint mass center j;

joint inertia tensor expressed at
joint mass center j;

generalized member stiffness matrix;
length of a member;
joint mass matrix expressed at joint J;

joint mass matrix eXpressed at joint
mass center j;

torsional moment;

fully plastic torsional moment;

bending moments;

fully plastic bending moments;

mass of branch B;

mass of joint J;

axial stress resultant;

subscript to denote negative end of a member;

position vector of member end N expressed in
global coordinate system;

subscript to denote section x = (l-B)L
of a member;
subscript to denote section x = (l-28)L

of a member;
subscript denoting positive end of a member;

position vector of member end P expressed in
global coordinate system;

external joint force (vector);

internal joint force (vector);

viii

N
No

force (vector) on joint J;

force (vector) on joint mass center j;

subscript to denote section x BL of a member;
subscript to denote section x = 28L of a member;
superscript denoting "plastic part";

rotation matrix (6x6);

rotation matrix (3x3);

elastic member stiffness matrix;

geometric member stiffness matrix;

submatrix of 86;

force transformation matrix defined by
Equation (l-l);

time;

superscript denoting matrix transposition;
member displacement (vector);

linear member displacements;

yield surface normal in dimensional force
space

member distortion (vector);

position vector of branch mass center b
expressed in global coordinate system;

position vector of joint mass center j
expressed in global coordinate system;

position vector of joint mass center j
expressed in local coordinate system;

joint displacement, velocity, and
acceleration (vectors), respectively;

acceleration (vector) of joint J;

ix

acceleration (vector) of joint mass
center j;

foundation acceleration (vector);

joint (global) coordinates;

member (local) coordinates;
responsiveness matrix for joint J;
dimensionless parameter modifying
moments and products of inertia of

a joint rigid body;

dimensionless parameter giving plastic-
portion length of a yielding member if
multiplied by 2L;

dimensionless parameter modifying
direction of yield surface normal;

prefix denoting "change" or "increment";
axial shortenings of a member divided by L;

uniformly distributed mass per unit length
of a member;

density of material;

yield function;

elliptic yield function

parabolic yield function

axial force parameters;

rotational member displacements; and

prefix denoting "gradient."

 

orcb;.
.
finn‘.‘
I

U“‘I¢.

enoug,

 

 

 

“WP:

takih

 

CHAPTER I

INTRODUCTION

This chapter presents the objective of the present
work, previous related studies, and assumptions and limita-
tions of the analysis develOped. It also gives an outline
of the investigation carried out and certain general defi-

nitions needed in the subsequent analysis.

1.1 Objective

 

Two types of nonlinearities occur in structural
problems. The first type may be referred to as "geometric
nonlinearities." They occur when deflections are large
enough to cause significant changes in the geometry of the
structure. In this case the equations of equilibrium must
be formulated for the deformed configuration. The second
type may be referred to as "material nonlinearities" which
include any deviation from linear elasticity, such as non-
linearly elastic, or plastic, or viscoelastic behavior of
the structural material.

The objective of this thesis is to develop a
numerical method for the dynamic analysis of space frames,

taking the effects of geometric and certain material

 

wor}

des:

1.2

haVe
the

Pape
recei
bys

and

the

 

ment

list

 

nonlinearities into consideration. In order to accomplish
the objective, a matrix formulation of the problem is de-
rived. Furthermore, a computer program is prepared for
the implementation of the analysis. Finally, numerical
results of certain problems are obtained in order to demon-
strate the validity and practicality of the method.

It seems hardly necessary to point out that the
problem under consideration is a broad, and hence, diffi-
cult one. However, as it will become apparent later, the
present work has had the advantage of using several earlier
works as "building blocks." These earlier works will be

described briefly in the following section.

1.2 Previous Studies

Recent advances in the field of computer technology
have provided the necessary tools for the development of
the analysis of geometrically nonlinear frames. Many
papers concerned with such developments have appeared in
recent years. Among these perhaps the pioneering works
by Saafan (12),*Idvesley (7), Argyris (l), and Johnson
and Brotten (6) deserve special attention.

In a recent paper, Jenning (5) has incorporated
the effects of change of geometry into certain displace-

ment transformation matrices for members of plane frames.

 

*Numbers in parentheses refer to entries in the
list of references.

He considered the axial shortening due to member inclina-
tion as being the only important nonlinear term in cases
where displacements are not exceedingly large. However,
he mentioned that accurate transformation matrices can be
adopted in cases in which very large deformations are to
be dealt with. Based on Jenning's formulation, Iverson
(4) has derived stiffness coefficients for space-frame
members. He applied his geometrically nonlinear formula-
tion to elastic frames under dynamic loads.

Connor, Logcher, and Chan (2) have also derived
a geometrically nonlinear formulation for the three dimen-
sional case. Their derivation is restricted to the small
rotation case, i.e., a case in which the squares of rota-
tion angles are negligible with respect to unity. Zarghamee
and Shaw (15) have given a similar formulation independently.
They presented a more comprehensive expression for axial
force than the one given by Connor, Logcher, and Chan.

The area of inelastic behavior of frames under
static and/or dynamic loads has been of great interest
in the past decade. Almost all published works on this
subject have been partially based on the plastic potential
theory. In a recent report, Morris and Fenves (9) have
studied the inelastic behavior of space frames under static
loads. They derived the incremental force-deformation re-
lations for a member having any number of plastic hinges.

In their report, they also derived the approximate yield

surface equations for certaincommonly used cross-sections,
considering the interaction between bending moments, tor-
sional moment, and axial force.

Nigam (11) has recently presented an elastoplastic
formulation to study the response of dynamically loaded
space frames. He derived the incremental force-displacement
relations for a member whose one or two end cross-sections
are yielding. His equations are basically the same as the
ones derived by Morris and Fenves, although they apparently
seem different. He then extended his formulation to the
dynamic case. But, in this extension, he did not discuss
such important questions as: how to handle the mass; how
to prevent the force vector acting on a yielding cross-
section from proceeding beyond the yield surface; and how
to solve the equations of motion.

Wen (13) has also studied the elastOplastic be-
havior of space frames under dynamic loads. His formula-
tion is based on dividing a yielding member into an
elastic member with continuous flexibility and an elasto-

plastic member with lumped flexibility.

1.3 Assumptions and Limitations

The assumptions and limitations employed in this
study are divided into the following five categories:
(1) Material
(a) The material is assumed to be linearly elastic——

perfectly plastic.

(b) All stress-strain characteristics of the
material are assumed to be time-independent.

(c) Strain-hardening effects of the material are
neglected.

(2) Cross sections

(a) All cross-sections are assumed to have two
axes of symmetry.

(b) Yielding is assumed to occur at individual
cross-sections with no "spread length."

(c) A shape factor of 1.0 is assumed for all
cross-sections. That is, cross-sections are
assumed to make an abrupt transition from a
completely elastic state to a state where
unrestricted plastic flow can occur.

(d) It is assumed that cross-sections are free
to warp under torsional loads, i.e., the case
of "pure" torsion prevails.

(3) Members

(a) All members are assumed to be prismatic.

(b) The usual engineering theory of bending is
assumed to be applicable.

(c) Torsion-flexure coupling effects are neglected.
(d) Members may have uniformly distributed gravity
loads only. Concentrated loads at points

other than the ends of members are not per-

mitted.

....- V
_- I.

f.“

"H
,.
'1
J?

_|

 

 

 

1.4

 

(4) Joints
(a Joint loads, which are necessarily concentrated,
may be statical and/or dynamical.
(5)FTmM3
RN Changes in the geometry of the frame are re-
stricted to the "small rotation" case, i.e.,
a case in which the squares of rotation angles

are negligible in comparison with unity.

1.4 Outline of Present Study

 

In this study, two yield condition equations for
a member cross-section, as used in Reference 13, are first
specified. Then a series of incremental force-displacement
relations for a member, taking the effects of geometric
onlinearities into account, are derived. Firstly, these
elations are summarized for a linearly elastic member,
; fornuLLated by Connor, Logcher, and Chan (2). Secondly,
:e case of a member whose end cross-sections are yielding
studied. Finally, the incremental force-distortion
lations for a discrete model, taken from Reference 13,
2 derived.
TTue dynamic analysis begins with a mass lumping
:edure (see References 4 and 13) that accounts for
try inertia. Then, the motion equation for a typical
(ungrounded) joint is derived. In the dynamic solu-

, the criterion by which a member cross-section is

 

 

 

y—y-w.

j...

 

 

rule;
next
secti
face
dyna:
of m“
vide.

gene.

men:

1.5

 

 

 

whos
Coun
0f

ruled to be yielding is described first. It is shown
next that the force vector acting on a yielding cross-

section is prevented from proceeding beyond the yield sur-

face. Finally, the transient response {of the frame to

dynamic loading is obtained by integrating the equations

of motion numerically. Dynamic loading is either pro-

vided by concentrated loads applied to free joints, or

generated by ground motion due to an earthquake.

Three numerical examples taken from Reference 13

are considered: a cantilever beam subjected to a pulse

type of loading; a six-member space frame subjected to a

step-function type of loading; and a two-story two-bay
building frame subjected to the 1940 El Centro earthquake.

These examples, on the one hand, illustrate the applications

of the numerical method developed. On the other hand, they

will provide a basis for comparing the numerical results

of the present method with those reported in the above-
mentioned reference.

1. 5 General Definitions

Figure 1-1 illustrates a threebdimensional frame
whose joints and members are all numbered. The joint

count is separate from the member count. For convenience

of computer programming, free joints are always numbered
first.

Each member is arbitrarily assigned an orientation

by specifying one of its ends as the positive end and the

other as the negative end. Consider a member framed be-

tween two joints. It is then said that the member is

positively incident to the joint at its positive end and

negatively incident to the joint at its negative end. For

example, in Figure 1-1, the member 2 is positively incident

on the joint 1 and negatively incident on the joint 3.

Two classes of right-handed Cartesian coordinate

systems are used. The origin of each coordinate system

is either clear from context or specifically pointed out.

The first class consists of a single joint (global) coor-

dinate system arbitrarily chosen. The second class con—

sists of all member (local) coordinate systems chosen in

the following manner. The first axis of each member coor-

dinate system coincides with the undeformed centroidal

axis, and is directed from the positive end to the negative

end of the member. The other two axes are coincident with

the principal axes of inertia of a generic cross-section

of the member .

In the next two chapters, force and displacement

transformation matrices will be needed. These are given

below for the sake of completeness. Referring to Figure

1-2, let olxlylzl and ozxzyzz2 denote two coordinate systems

whose axes are parallel and oriented in the same fashion.

The force transformation matrix from the first system to

v Lani . .1

 

 

t t
.J
(1)

RE fe;
Coorc
VI

02):

SYSte

Whery

 

the second is then given by

(1-1)
H I .

where I is the third order identity matrix, and

F—

 

 

0 Z -Y
H = -Z X (1-2)
a

iJi vfluich X, Y, and Z are the coordinates of the center 02
veitfi respect to the system lelylzl. The corresponding
displacement tranformation matrix is given by

-1 t
D = (T ) _
1 2 0201 (1 3)

Referring again to Figure 1—2, let Ozxéyézé denote a third
coordinate system rotated with respect to the system
02x2y222. The force transformation matrix from the second

system to the third is given by

 

r 0
R = (1—4)
0 r
where
12g5(xé,x2) cos(xé,y2) cos(xé,22)
_ I l _
r — cos(y2,x2) cos(yé,y2) cos(y2,22) (l 5)
I O O
L:os(zz,x2) cos(22,y2) cos(zz,22)

 

_J

10

in which the parenthesized quantities denote the angles
between the axes indicated. The matrices R and r will

be referred to as the "rotation" matrices. In this case,
the displacement transformation matrix is also given by
the matrix R. For example, the incremental member end
displacement AU and the incremental displacement AX of

the joint, connected to the member end being considered,

are related by
AU = RAX (1-6)

where AU and AX are expressed, respectively, in the member
and joint coordinate systems, and R is the rotation matrix

from the second system to the first.

 

w-“ -. v I

 

N

1

 

 

 

CHAPTER II
FORCE-DEFORMATION RELATIONS FOR A MEMBER

In this chapter, two yield condition equations for
a member cross-section, to be employed for obtaining nu-
merical solutions, are first specified. Then a series of
incremental force-displacement relations for a member,
considering the effects of geometric nonlinearities, are
derived. Firstly, these relations are summarized for a
linearly elastic member. Secondly, the case of a member
whose end cross-sections are yielding is studied. Finally,
the incremental force-distortion relations for a discrete

model, as used in Reference 13, are derived.

2.1 Yield Condition Equation

In frame analysis, it is convenient (perhaps even
necessary) that the yield condition equation for a member
cross-section be formulated in terms of stress resultants.
These resultants will be hereafter referred to as "force
components" for short. EXpressed in these terms, the
yield condition equation defines the combination of the
force components necessary to initiate yielding at a

cross-section. Since a shape factor of 1.0 is assumed

for all cross-sections, the initiation of yielding will

11

12

coincide with the start of unrestricted plastic flow.
Thus, the yield condition equation represents the relation
among the force components for the initiation of unre—
stricted plastic flow. Finally, for an elastic-perfectly
plastic material as assumed in this study, the yield con-
dition equation remains the same as yielding progresses.
In the case of three-dimensional frames, the yield
condition equation for a cross-section may be written

symbolically as
¢(FX,MX,My,MZ) = 1 (2-1)

where 4 is the so-called yield function, Fx is the axial
force, Mx is the torsional moment, and My and M2 are the
bending moments. Equation (2-1) represents a hypersurface
in the four-dimensional force space spanned by the Cartesian
coordinates FX, Mx’ My' and M2. This hypersurface is called
the yield surface.

The derivation of the exact yield condition (or
yield surface) equation for a given cross-section is gen-
erally quite difficult. Numerical solutions are available
only for a square cross-seciton. Morris and Fenves (10)
have derived approximate lower bound yield surface equa-
tions for commonly used cross—sections, using a procedure
suggested by Hodge (3) along with simplifying assumptions
on the neutral axis position. Even those equations appear
rather unwieldy for general application to frames with

many members.

 

the fur;

 

13

In this study, for the sake of simplicity, two

particular forms of yield functions taken from Reference

13 are employed. These functions do not correspond exactly
to any commonly used cross-section. However, considering
the assumptions made and uncertainties involved in applying
the theory of plasticity to this problem, they may be re-
garded as reasonable approximations for obtaining numerical
solutions.

The first function is given by

 

(u'

 

[FX 1 Mx 2
¢p = EFr—' + (M——)
I XP XP
’M _ Mz
+ :M-L: + _M = l (2_2)
i YP, ' 2P

 

in which Fxp' M , M , and sz are the fully plastic

YP

force components corresponding to Fx, Mx’ My' and M2,
respectively. The function ¢p will be referred to as the
"parabolic" yield function. The second function is given

by

F 2 M 2
—_’_‘_ .3‘__
ee-(P>+(M>
XP xP

fr. 2 M2 2
+ (M l + (571—) = 1 (2-3)
YP 2P

the function @e will be referred to as the "elliptic"

yield function.

p4

‘ a
Irv-n ' 4.?
“yob-la¢ i:

qs- .‘
”H. n
vaoor- \-

I
2 v...
to ole“,

~e .' 1
1:.3V a:

vb. ‘

u
'0 "

. "A
.'.a"~~ln

 

14

It is apparent that the parabolic yield surface
is enclosed by the elliptic one. Therefore, for engineering
analysis, it is a more conservative one to use as it would
normally indicate an occurrence of yielding before the
elliptic yield function would do so.
2.2 Incremental Force-Displace-
ment Relations for a

Geometrically Nonlinear
Member

 

 

 

The incremental force-displacement relations for
a typical member can be derived by considering the statics,
member geometry, and stress-strain characteristics of the
material. Connor, Logcher, and Chan (2) derive these
relations for an elastic member, considering the effects
of geometric nonlinearities. This derivation will be
outlined briefly in the next subsection for the sake of
completeness.

When a member end cross-section yields, an un-
known vector, namely, the plastic end displacement incre-
ment, is introduced in the incremental force-displacement
relations. This unknown is determined, up to a scalar
parameter, i.e., the "flow constant," by applying the
plastic potential theory. (This theory is described in
several textbooks among which the one by Malvern (8) is
perhaps the most recent.) The flow constant, and con-

sequently the force-displacement relations, can be

:.. rm a;
i:- .ou.‘l\.

 

 

2.2.1

I r”

Of a tv

' (,1

. p ’
nae,” fl
‘Bov‘ane“

' a
Cfl‘na. q
‘V‘ao'.“

 

bible.

 

 

15

determined by applying the condition that the force incre-
ment--to be theoretically correct, of infinitely small
magnitude--acting on a yielding cross-section is tangent

to the yield surface.

2.2.1 Elastic Member

 

Figure 2-1 shows the initial and deformed positions
of a typical member PN. The x axis coincides with the un-
deformed centroidal direction and the y and z axes are
coincident with the two principal axes of inertia. These
axes constitute the local coordinate system for the member.
M

F

The internal forces (F M , M2) and displace-

y’ X' y
wz) are referred to in this

2!

ments u u u
( x! Y, 2' (0X, (Dy!

coordinate system, and so are the end forces FP and PN
F

with components FP etc., and end displacements UP'

Nx’

and UN with components U

x!

U etc.

Px' Nx'
Figure 2-2 shows the initial and deformed positions
of a differential element AB. The centroidal point B is
displaced to the point B‘. Based on the assumption of
small rotations (w; << 1, w: << 1), the longitudinal strain

e at the point B', in the direction of the tangent to the

deformed centroidal axis, is expressed by

du

_ x 2 2 _
e — di— + 0.5 my + 0.5 oz (2 4)

'6" 3'0

H...»

 

 

with
du
w _..__z.
y dx
du
wz = dx
where u , u , and u
X Y

the point B.

The second and third terms in Equation

16

(2-5a)

(2-5b)

z are the translational displacements of

(2-4)

. -q, -_, ,-~

introduce the effects of geometric nonlinearities in the

”A in-..“ :4

.- —~.-.-' urn - "

present formulation. If these terms are neglected the

.. ~_r-

formulation

“#0” -~- -

Mo-—-_.-M~.

.. c,.

returns to the geometrically linear case.

With reference to Figure 2-2, the force (stress

resultant) tangent to the deformed centroidal axis is

denoted
used to
meaning
obvious

related

by

by N.
denote the negative
in every case where
The

from context.)

to the strain e and

(Note that the same symbol N has also been

end of the member. But the
the symbol N appears should be
force N and the moment Mz are

displacement uy, respectively,

(2-6)

(2-7)

in Which E is the modulus of elasticity, AX is the cross—

seCtional area, and I2 is the cross-sectional moment of

inertia

about the z axis.

Note that Equation (2-7)

is

 

based or

neglects

present
he Obti
The de:

straig'r.

Equatio

along t,

 

Equati

member.

in whic

17

based on the usual engineering theory of bending which
neglects the effects of shearing deformations.

The equilibrium equations consistent with the
present geometric approximation (w; << 1” w: << 1) can
be obtained by applying the principle of virtual work.

The details are omitted since the procedure is quite

straightforward. The final equations are then listed as

§§.= o (2-8)
2
d M

Equation (2-8) states that the axial force N is constant
along the entire length of the member. Consequently,
Equation (2-6) can be integrated over the length of the

member. This leads to

EAx

N = - -L— (qu-uNx) + EAx (6y+62) (2-10)

in which L is the length of-the member and

..1__ IL 2
5y 2L 0 “zdx
(2-11)
L
= 1. I 2
52 ‘21. o “’d"

+‘—
5.!
e

 

cte
-Vov

.\

:(‘Vn

 

I

Q..;Cv

‘V

. .
n‘qfiv‘"

«“

J-f‘

he».

a3
. \L

,-

one Q'

as S u‘T'lL“:

3a
a

D;

18

The governing equation for the displacement uy is
obtained by substituting Equations (2—5b) and (2—7) into

Equation (2-9). Thus

 

d4u ¢2 d2u
+ 2y 22: 0 (2—12)
dx L dx
where
0.5
2
NL
¢y - if” (2-13)

axial force parameter

Note that the parameter ¢y is a function of the axial
force N which in turn depends on the displacement uy
Rigorous solution of the differential equation (2-12) for
the displacement uy would require using iterative schemes.
To avoid that, the latter equation should be solved under
the assumption that the parameter ¢Y (or axial force N)
is known. The solution should then be introduced into
Equations (2-7) and (2—10). The discrepancy between the
assumed axial force N in solving Equation (2-12) and the
one obtained in Equation (2-10) is of course a source of
inaccuracy in the procedure. This point would be brought
up again in later sections.

It is evident that Equations (2-7), (2—9), and
(2-12) are associated with bending in the x-y plane.

Similar equations can readily be written for bending in

 

twisti:

mined i'

in whi
metric.
linear
matr13.
effeCt
haVe 5
it has
releas
II) t1:

menta;

 

in in:

19

the x-z plane. The force-displacement relations are ob—
tained by assembling the equations associated with flexure
in the x-y plane, flexure in the x-z plane, and finally,
twisting about the x axis. The end forces are then deter-
mined in terms of the end displacements from the latter
equations by enforcing the prOper boundary conditions.
The lengthy details of this procedure are omitted. Only
the relevant equations in matrix form will be given below.
The end force-end displacement relations are now

expressed as

P PP P PN N G
t (2-14)
FN = SPN UP + SNN UN + FG
' ' ll ° ll '
in which SPP' SPN' and SNN are the elastic stiffness

matrices, and FG is a column matrix containing the non-
linear terms due to the rotations wy and wz. These
matrices are all expanded in Appendix I. Note that the
effects of torsion-flexure coupling and warping restraint
have been neglected in Equations (2-14). In addition,
it has been assumed that the member does not have any
releases. However, it will be shown later (see Appendix
II) that end releases may be considered in the incre—
mental form of the latter equations.

Later, the end force-end displacement relations

in incremental form will be needed. To derive these, it

20

S and S

is assumed that the stiffness matrices S PN’ NN

P P '

are constant within the increments (of forces and dis—
placements). This assumption is a reasonable one since

the elements of S S and S vary slowly with oi

PN' NN

(i = y,z) when ¢i is not close to Zn. (¢i = 2n corresponds

PP'

to N = the Euler load.) The incremental relations are

thus expressed as

AFP = SPP AUP + SPN AUN - AFG
AFN = SPN AUP + SNN AUN + AFG
where the increment AFG can be expanded as
AFG = SG (AUP - AUN) (2-16)

The matrix SG may be interpreted as the "geometric" stiff-

ness matrix. The approximate form of this matrix used

in the present analysis is developed in Appendix I.
Introducing Equation (2-16) into Equations (2—15),

the approximate incremental relations are finally written

as

AFP = KPP AUP + KPN AUN
t (2-17)
AFN = KPN AUP + KNN AUN

in whi':

 

 

 

21

in which

KPP = SPP + SG
KPN = SPN - sG (2-18)
KNN = SNN + so

It follows from the symmetry of the stiffness

matrices S S and S that the generalized stiffness

PP' NN' G

matrices K and KN

PP are also symmetric.

N
If the member happens to have end releases, its
incremental force-displacement relations must be modified
accordingly to account for such releases. This modifica—
tion can be carried out easily in a fashion quite similar
to that presented in Reference 9. However, for the sake

of completeness, an outline of the latter modification

is given in Appendix II.

2.2.2 Elastoplastic Member

 

The assumption that the member is elastic is now
dropped. It is instead assumed that both end cross-
sections of the member are yielding. The end displace-
ment increment AU can be decomposed into an elastic part
and a plastic part. These two parts will be denoted by

the superscripts e and p, respectively. Thus

in whiz

-‘l

functiJ
at the
A gee:
in Fig
COndit

Seen 5

 

 

 

 

 

22

AU = AU

8 P
P P + AUP

(2-19)

AU = AU

9 P
N N + AUN

The above decomposition is a basic assumption of the theory
of plasticity.
Based on the plastic potential theory, the flow

law adapted to the present problem is expressed by

P = = A
AUP bP V¢ (FP) vaP
(2-20)
p _ _
AUN - bN V¢ (FN) — vaN

in which b is the flow constant (a positive scalar parameter),
V is the so-called "del" (gradient) Operator, ¢ is the yield
function, and V is the outward normal to the yield surface
at the point where the end force F meets the yield surface.
A geometrical interpretation of the flow law is illustrated
in Figure 2-3 for a simpler case in which the yielding
condition depends on only two stress resultants. It is
seen from this figure that the plastic end displacement
increment AUp is normal to the yield surface. Thus, the
direction of the increment AUp is determined by the plastic
potential theory. However, the magnitude of that incre-
ment still remains to be determined.

During yielding, the end force vector F must stay

on the yield surface represented by the equation<b(F) = l.

(The latter equation remains fixed since the effects of

or,

:hn L '

“1“3»1C
mignlt;
M"

Llese e

tngent

Per-feet

by the

fOre, E

 

23

work hardening are neglected.) This is specified analyti-

cally by
A¢ (FP) = 0
(2-21)
A4 (FN) - 0
These equations can be rewritten as
[VcMF Ht AF = o
P P
t (2-22)
[V¢(FN)] AFN = 0
or,
t _
VP AFP — 0
t (2-23)
VN AFN = 0

Equations (2-23) will be used shortly to determine the
magnitude of the incremental plastic end displacement AUP.
These equations state that the end force increment AF is
tangent to the yield surface.

Since the material is assumed to be elastic--
perfectly plastic, the end force increment AF is governed
by the elastic and displacement increment AUe only. There—

fore, Equations (2-17) can be written as

_ e e

APP ‘ KPP AUP + KPN AUN
(2-24)

_ t e e

AFN ‘ KPN AUP + KNN AUN

24

Considering Equations (2—19) and (2-20), the above

equations can be rewritten as

AFP = KPP AUP + KPN AUN
‘ KPP VP bp ' KPN VN bN
t (2-25)
AFN = KPN AUP + KNN AUN
- Kt v b - K v b
PN P P NN N N
At this point, it is convenient to define two
6x2 matrices GP and GN by
I
GP ‘ [KPP VP 1 KPN VN]
t | (2-26)
GN z [KPN VP I KNN VN]

If now the flow constants bP and bN are assembled into a

column matrix B, Equations (2-25) can be expressed as

AFP = KPP AUP + KPN AUN - GP B
(2-27)
AFN = KPN AUP + KNN AUN - GN B
Substitution of these equations into Equations (2—23)
furnishes
t t t _
VP KPP AUP + VP KPN AUN VP GP B — O
t t (2-28)
v t AU _ t =
N KPN P + VN KNN AUN VN GN B 0

These a

estat;

eh u ..
Ziernt'.

in w'r~

A-

:P“ae-
u“ -4
.

beco:,

 

 

 

25

These equations are now combined into the following single

equation:

"t - t - r-t -
VP KPP FY? KPN VP GP

AU + AU - B = 0 (2—29)
vt Kt P vt K N vt G
LN P12 -N Na -N 1!.

 

 

 

 

 

 

Solution of Equation (2-29) for the (flow constant)

matrix B yields

t t

B = E GP AUP + E GN AUN (2-30)
in which
r- -1 ‘1
t
VP GP
N N
L. _J

 

 

Equation (2—30) is finally substituted back into Equations

(2-27) whereby the incremental force-displacement relations

become
AF = (K - G E Gt) AU
P PP P P P
+ (K - G E Gt) AU
PN P N N
(2—32)
_ _ t t
AFN — (KPN GP E GN) AUP
+ (K - G E Gt) AU
NN N N N

ab; rv
seVV‘AV-

 

M: H»-
auu we»

Spond

Elast

 

26

The reduction to the case where only one end cross-
section of the member is yielding is quite straightforward

and thus omitted.

2.2.3 Elastic Return

 

When a cross-section reaches its yielding condition,
a "plastic hinge" is said to form there (following the
terminology in the simple plastic theory of structures).
It frequently happens that a plastic hinge "unloads" and
the corresponding cross-section becomes elastic again.
This phenomenon is referred to as an "elastic return.“ It
occurs whenever there is a reversal in the direction of
the incremental plastic displacement at a yield hinge.

It thus follows from Equations (2-20) that the elastic
return at a plastic hinge is signalled by a negative flow
constant.

To test for elastic returns, the column matrix B
for each elastoplastic member is first calculated. The
sign of each of the elements (one or two) of the matrix
B is tested. Any yield hinge at which an elastic return
is detected is then removed; in other words, the corre-
Sponding member end cross-section is considered to be

elastic.

2.3

Itfinn‘l_L1|*4

   
  
  
  
  

menbe
direc1
hwevé
cally
this 1

disto

 

isd

to t1-

 

27

2.3 Incremental Force-
Distortion RéIatIEns for
a Discrete_Model ofia 530-
metricaIly Linear MemBer

 

It is evident that for a "geometrically linear"
member the force-deformation relations can be obtained
directly from the previous results. These relations are,
however, rederived in the present section for a geometri-
cally linear discrete model taken from Reference 13. To
this end, the derivation will have to be made in terms of
distortions rather than displacements.

Figure 2-4 shows a model PN which is to replace a
typical geometrically linear member PN whose P1 and N1
sections are yielding. It is assumed that each of the
segments PPl, Ple, NNl' and NlN2 (of the member) are
infinitely rigid and of length BL where 3 is a nondimen-
sional parameter to be chosen between zero and some
fraction, say 1/8. Furthermore, it is assumed that the

flexibility of the portion PP is lumped at its midsection

2

P The same assumption is also made for the portion NNZ'

1.
The remainder of the model, namely, the portion PZNZ is
assumed to be continuously elastic.

The distortion of the model referred to the end P
is defined to be the displacement of the end P relative

to the end N. In incremental form, it is given by

= AU - D AU (2-33)

AW P PN N

PN

 

 

in W111

iron ‘

 

distc‘

The p

If th

col

in wt

The e

Sunny

 

28

in which DPN denotes the displacement transformation matrix

from the local coordinate system through N to the same
system through P.
Analogous to Equations (2-19), the incremental

distortion AW can be decomposed into

PN

= e P _
AWPN AWPN + AWPN (2 34)

The plastic distortion increment AWPN can be expressed as

P - _
AwPN — DPP1 VII bpl + D v b (2 35)

If the flow constants bP and b are assembled into a
l 1

column matrix B, Equations (2-35) can be written as

N

P 3 _
AWPN GB (2 36)
in which
G = [D v I D v ] (2-37)

The elastic distortion increment AWe can be given as a

PN
summation of

6 e e
AW = D AW + D AW
PN PPl PlPl PP2 PZNZ
+ D Awe
PN N N (2-38)

1 1 l

 

in wk;

gortL:

 

33m ‘
17 ‘ a
515.6.

matrip

Where

 

 

29

The flexibility of the portion PP2 is lumped at its mid—

section Pl; therefore, the elastic distortion increment

Aw: P is determined by
l l
e
AW = AP AP = AP T AF (2-39)
P1P1 1Pl P1 lPl P1P P

in which AP P is the lumped flexibility matrix for the
l 1

portion PP2, and TP P is the force transformation matrix
1

from the local coordinate system through P to the same
system through P1. The main diagonal elements of the

matrix AP P (which is diagonal) are given by
l

1
_ 28L _ 28L
C1 ‘ EA C2 ‘ GA
X Y
_ 28L _ 28L _
C3 ‘ 3A2 C4 - SJ (2 40)
_ 28L _ 28L
C5 ‘ E C6 ‘ EI
Y z

in which G is the modulus of shear, Ay and A2 are the
effective areas of shear, and J is the constant of torsion.

The increment AW: N is similarly determined by

e
AW = AF = T AF (2—41)
NlNl ANlNl Nl “NlNl NlP P
where
AN = A (2—42)
1N1 P1P1

30

Since the segment P2N2 is elastic, the superscript e can

be omitted from the increment AWe This increment is

P2N2

readily given by

AW = aw
P2N2 P2N2

= AF
APzpz P2

T AF (2—43)
AP2P2 P2P P

where AP P is the ordinary flexibility matrix for the
2 2

continuously elastic segment P referred to the section

2N2

P Note that the effects of shearing deformations may be

2.
included in the matrix AP P
2 2'

Equations (2-39), (2—41), and (2—43) are substituted

back into Equation (2-38). This leads to

e

AWPN = APP AFP (2-44)
or,
_-le _
AFP — APP AWPN (2 45)

(2-46)

f\
”V

 

31

Considering Equations (2-34) and (2—36), Equation (2-45) is

cast into

_ -l _
AFP — APP (AwPN GB)

Analogous to Equations (2-23), the

force increments AFP
1

surface is expressed by

and AF

(2-47)

condition that the

are tangent to the yield

1

- GB)

These two equations can be combined into

Gt Ag; (AWPN - GB)
from which
B EGt Ag; AWPN
where
1

E = (Gt Ag: G)-

0

0
(2-48)
0
(2-49)
(2-50)
(2-51)

Equation (2-50) is substituted back into Equation (2-47)

to furnish

AFP =

-l -1
(APP ’ APP G E G

t

APP

) AW (2-52a)

PN

32

The incremental force-distortion relations are finally

completed by having

AFN = --TNP AFP . (2-52b)

CHAPTER III

DYNAMIC ANALYSIS AND SOLUTION

In this chapter, the joint mass matrix is first
constructed. The equation of motion for a given joint is
derived next.' Finally, the numerical procedure employed

to determine the dynamic response is presented.

3.1 Dynamic Analysis

 

The dynamic analysis begins with the formulation
of the equation of motion for a typical free (ungrounded)
joint in the global coordinate system. This requires that
a mass matrix be constructed for the joint. Such a con-
struction may be carried out by lumping, at the joint,
the "contributory" masses from all the incident members.
The mass lumping could be done in a reasonable manner as
long as there is no coupling between the various joint
mass matrices. However, the more realistic the lumping
procedure is, the more accurate the numerical results will
be.

In this study, the joint mass matrix is constructed
in a manner described by Iverson (4). The rotary inertia
is taken into account in this construction. A nondimen-
sional parameter, as introduced in Reference 13 to control

the moments and products of inertia, is also included.

33

dairiekflj

34

3.1.1 Joint Mass Matrix

 

To formulate the joint mass matrix, it is convenient
to envisage a rigid body associated with the joint. To this
end, each member with its length multiplied by a will be
called a "branch," where a is a dimensionless parameter, to
be chosen between zero and 1/2. The joint rigid body is
then defined to be the collection of all the branches incident
to the joint being considered. To keep the total mass of
the rigid body independent of a, the density p of the
material is divided by 0. Furthermore, it is recalled
that each member may have a uniformly distributed gravity
load with a mass u per unit length. Then each u is also
divided by a. It follows that a is introduced to modify
the moments and products of inertia of the joint rigid
body while leaving the rigid body mass the same.

The mass mB of a typical branch B is presented

readily by
mB = 0.5 L (Ax p+u) (3-1)

The mass m of a generic joint J is obtained by summing the

J
masses of all the incident branches. Thus

mJ = E’mB (3-2)

It is seen that a does not appear in Equations (3-1) and

(3-2).

35

Referring to Figure 3-1, the position vector

 

(Xj, Yj’ Zj) = Wj of the joint mass center j in the joint
J coordinate system is determined by
Em X
X. = mB b
3 J
Zm Y
y. ==——35133 (3-3)
3 J
2m Z
Z. = mB b
3 J

in which (Xb, Yb' Zb) = W5 is the position vector of the
branch mass center b in the same coordinate system, and
each summation is to be taken for all the branches incident

to the joint J. The vector W5 is given by

wb = e 0.25 o (P - F) (3-4)

in which P and N are, respectively, the position vectors
of the positive and negative ends of the member B in the
global coordinate system. The positive or negative sign
in Equation (3-4) is to be chosen according as the member
B is positively or negatively incident to the joint J.
The branch inertia tensor Ib with respect to the
local coordinate axes through b is determined next. Note
that these axes are coincident with the principal axes

of inertia for the branch B. Thus

1b = 0 C2 0 (3-5)

_. C111

 

 

36
where

C = I”p (y2 + 22) dx dy dz

F‘

_ 0
-7 (1y + 12) (3-6a)

C = 1110 (z2 + d2 xz) dx dy dz

+11.) (0.2 X2) dX

2 3

a L

C = 111p (y2 + a2 x2) dx dy dz

+ Iu(a2 x2) dx

2 3

a L 0L _
—-—-96 (p Ax + p) + I (3 6C)

7 y

The tensor Ib is now transferred to the joint

mass center j by using the theorem of parallel axes. This

 

 

leads to
I] = Ib
ryz + z2 -x y -x z-fl
j j j j J J
2 2
+ m z. + x. - . 2. (3-7)
B J J y] 3
X2 + 2
symmetric j yjj
L—
where (xj, yj, zj) = Wj is the position vector of the joint

mass center j in the local coordinate system through b.

37

This vector is determined by

W. = r (W. - W

in which r is the rotation matrix (3x3) from the joint J
coordinate system to the member B coordinate system. The
definition of the matrix r is given by Equation (1-5).
Each branch inertia tensor Ij at the joint mass
center j is now rotated to the global coordinate system
through j. They are then summed to obtain the joint

inertia tensor Jj at the joint mass center j. Thus
J. = Z r I. r (3—9)

Finally, the mass matrix for the joint J with
respect to the global coordinate system through j is

assembled as

M. = (3-10)

where I is the third order identity matrix.

3.1.2 Equation of Motion for
a Joint Rigid Body

 

 

Considering again the joint J and its mass center
j in Figure 3-1, the equation of motion for the corre-
sponding rigid body in the global coordinate system

through j is eXpressed as

Mj x. = F (3—11)

38

in which Xj is the acceleration (vector) of j and P3. is
the resultant of all non-inertial forces acting on j.

This equation can be rewritten as

-1 t .. _
Mj (TjJ) XJ - TjJ PJ (3 12)

in which TjJ is the force transformation matrix from the
global coordinate system through J to the same coordinate

system through j. Premultiplying Equation (3-12) by

T71, it follows that
jJ
MJ XJ = PJ (3-13)
in which
t
g -1 -l
The inverse of MJ may be referred to as the "re-
sponsiveness" matrix. This inverse is denoted by ZJ and
given by
L. _
-l -l
I/mJ - HJ. H --HJj
-1
Z = M : _ _ (3‘15)
J J Jj1H le
L J

 

 

where H, as defined by Equation (1-2), contains the

coordinates (xj'Yj' Zj) of the joint mass center j in

the global coordinate system through J (see Figure 3-1).
The force PJ acting on the joint J, in the global

coordinate system through J, may be written as

39
P = P + P (3-16)

where PI is the internal force (vector) acting on J and
I

P is the external load (vector). The force P

E in terms

II
of the member end forces in their corresponding member
coordinate systems, is given by

(P) (N)

PI=-ZRFP-ZRFN (3-17)

in which the first summation is to be taken for all the
members positively incident to the joint J, the second
summation is to be taken for all the members with negative
incidence, and R is the rotation matrix (6x6) from the
apprOpriate member coordinate system to the joint J co-
ordinate system.

Summarizing, the acceleration of the joint J is

given by

X = Z (PI + PE) (3-18)

in which all the quantities are in the global coordinate

system through J. The responsiveness matrix Z and the

J
force PI are determined by Equations (3-15) and (3—17),
respectively. The external load P is to be prescribed.

E

If the foundation is subjected to translational
motion (e.g., due to earthquakes or blasts), it can readily
be shown (see Reference 13) that the equation of motion

becomes

+ P ) — M (3—19)

40

where i3 is now the acceleration of the joint J relative
to the foundation and i0 is the acceleration of the founda-

tion.

3.2 Dynamic Solution

The numerical method used to obtain the dynamic
solution will be outlined in this section. First, the
criterion by which a member cross-section is ruled to be
yielding is described. It will be shown next that
the force vector acting on a yield hinge is prevented from
proceeding beyond the yield surface. Finally the steps
of the solution procedure are described and the choice of
time increment is discussed.

3.2.1 Insertion of a Plastic

Hinge at a Cross-SectIon

It is presumed that a given cross-section is
initially elastic. This presumption is expressed analyti-

cally by
¢(F) < 1 (3-20)

in which 4 is the yield function and F is the force
acting on the cross-section. As the frame is deformed,
the force F may become so large that although Inequality

(3-20) holds, yet

0(F + AF) > 1 (3-21)

41

where AF is the force increment (say, corresponding to a
small time increment) computed under the assumption that
the cross-section is elastic.

If the above conditions hold, the cross—section is
ruled to be yielding. A yield hinge is inserted at the
cross—section and the force increment AF is recomputed
accordingly, i.e., with the newly inserted plastic hinge
taken into consideration.

3.2.2 Force Containment at
a Plastic Hinge

 

 

Let it be assumed that for some cross-section, a
situation as represented by Inequalities (3-20) and (3—21)
prevails, i.e., the cross-section is yielding. With ref—

(1)

erence to Figure 3-2, the unit normal V is determined

by

(l) - V¢(F) (3—22)

V ‘ |v¢(P)r

where V is the gradient operator. This unit normal is
used to compute the force increment AF. If the force

F + AF does not reach the exterior of the yield surface,
i.e., ¢(F + AF) 3 1, then the force increment AF is ruled
to be acceptable and no further action is taken. On the
other hand, if ¢(F + AF) > 1, then "force containment"
becomes necessary since it is physically impossible for

the force F + AF to proceed beyond the yield surface.

42

 

To this end, the unit normal V(2) is determined by
V(2) _ V4>(F + AF)
|V¢(F + AFN (13—23)

The force increment AF is then recomputed by using an

"average" normal V as given by
v = (1 - y) v”) +‘y V(2) (3-24)

where y is a scalar parameter to be chosen between zero
and unity. The choice of a value for y is discussed next.
Starting from zero, y is advanced by a finite value
up to the point where ¢(F + AF) is no longer greatern than
one. Limited numerical experience on the problems pre-
sented in the next chapter indicated that numerical results
were quite insensitive to the value of y. Thus, y was taken
to be either zero or one; in other words, if the unit nor-

(2)

mal V‘l) was not acceptable, the unit normal V was.

3-2.3 Numerical Procedure
of Solution

 

Assume that the state of the frame is known at a
given time t1. This includes: the joint displacements,
velocities, and accelerations; the member end forces;
and the list of the existing plastic hinges. It is then
desired to determine the state of the frame at the time
t = t + At where At is a finite time increment. To

2 1

this end, the following procedure is employed:

43

(1) Compute the incremental displacement of each free

(2)

(3)

(4)

(5)

joint from some numerical integration formula such

as
AX = [X(tl) + 0.5 X(tl) At] At (3-25)

Compute the incremental displacement of each member
end by rotating the proper joint displacement incre-
ment from the global coordinate system to the appro-
priate local coordinate system.

If the solution is to be geometrically linear, com-
pute the incremental distortion of each member from
Equation (2—33).

If the solution is to be geometrically nonlinear,
compute the stiffness matrices for each member
(elastic or elastoplastic) from Equations (A—l—8),
(A-1-9), and (A-l-lO). In this computation, the
unknown axial force of a member at the time t2 is
approximated by the known axial force of that mem—
ber at the time t1. It will be shown, in the first
example of the next chapter, that elimination of
such an approximation would not affect the numerical
results significantly.

Compute the incremental member end forces from
Equations (2—52) or (2—32) depending on whether

the solution is to be geometrically linear or

nonlinear. Then, increment the member end forces

by the values just obtained.

(6)

(7)

(8)

(9)

(10)

(ll)

44

Compute the yield-function value for each member
end force if the corresponding cross-section is
elastic (not yielding), by using such equations as
Equations (2-2) or (2-3). If this value happens
to be greater than unity, insert a plastic hinge
at the cross-cross being considered and return to
Step (5). (See Subsection 3.2.1.)

Remove every plastic hinge at which an elastic
return is detected. At this point, it seems
logical to return to Step (5) in order to recom—
pute the member end forces if necessary. But, to
avoid the possibility of removing and inserting a
particular plastic hinge repeatedly within a single
time-step, the return to Step (5) is neglected.
Compute the internal force acting on each free
joint from Equations (3-17).

Compute the acceleration of each free joint from
Equations (3-18) or (3-19) depending on the
absense or presence of the foundation motion.

Use such a numerical integration formula as

+ 0.5 [x(tl) + X(t2)] At (3-26)

to determine the velocity of each free joint.
Increment the free-joint displacements by the

values obtained in Step (1).

45

The state of the frame is thus completely determined at

the time t2. The same procedure is repeated for advancing

from the time t2 to the time t2 + At, and so on.

3.2.4 Choice of Time Increment

 

It is well known that a numerical procedure, as
applied to a linearly elastic structure, is stable if the
time increment used is smaller than a certain fraction
(say, l/n) of the smallest period of natural vibration.

This remark has to be ignored in the following discussion
since the system being considered is not linearly elastic.
However, knowledge of the smallest period is useful to the
extent that it serves as a guide in choosing the time incre—
ment. In general, the problem of determining the smallest
period is quite time—consuming. This problem becomes more
complicated by the fact that the smallest period will change
whenever a plastic hinge is inserted or removed. It would
thus appear more appropriate that the smallest period be
roughly estimated rather than rigorously computed. For

this purpose, the following procedure suggested by Iverson
(4) is employed.

The total mass of each free joint and the largest
axial stiffness of all the members incident to this joint

are used to compute

 

free-joint mass
axial stiffness

 

T = 2n

46

The smallest value of T for all free joints is taken as an
estimate for the smallest period. This estimate is then
used as an initial try for the time increment. In subse-
quent tries, increases or decreases are made if necessary.
In general, the largest tolerable time increment (that
yields a stable solution) should be used. Such a time
increment is normally so small that the corresponding
numerical results would not be significantly different

from those using smaller time increments.

3.3 Computer Program

 

A general computer program is prepared to imple-
ment the formulation presented. The program is written
in FORTRAN IV for the use on the CDC 3600 digital computer
at Michigan State University. It is described and also

listed in Appendix III.

’ J

CHAPTER IV

APPLICATIONS

This chapter presents three numerical examples
taken from Reference 13. These examples, on one hand,
illustrate the applications of the analysis developed in
the preceding chapters. On the other hand, they provide
a basis for comparing the present numerical results with
those reported in the above-mentioned reference. The
first example is a cantilever beam subjected to a pulse
type of loading with a short duration. The second one is
a three-dimensional rigid frame with a triangular plan-
form subjected to a step-function type of loading. The
final example is a two-story two-bay building frame sub—
jected to the 1940 El Centro earthquake.

In the examples presented, forces are expressed
in kips, and moments in kip—inches. The material is
assumed to be structural steel having the following prop-
erties: density = 490 pounds per cubic inch, Young's
modulus = 30,000 ksi, shear modulus = 12,000 ksi, normal
yield stress = 33 ksi, and shear yield stress = 18 ksi.
Linear displacements are expressed in inches, and rota-

tional displacements in radians.

47

48

4.1 Cantilever Beam

 

Figure 4-1 shows a cantilever beam which has a
uniform cross-section of 12WF53. The web of the beam is
vertical (in the X-Y plane). The fully plastic stress
resultants are 504, 275, 275 kips, and 73, 997, 2706 kip-
inches, respectively. Each of these values represents
the corresponding carrying capacity of the beam cross-
section if the cross-section is subjected to that one
stress resultant only. The external loads, in addition
to the weight of the beam, are given in the figure. Note
that the dynamic disturbance is supplied by a pulse type
of loading with a duration of 0.01 seconds. The graphs
to be presented in this example correspond to the elliptic
yield function, a value of a = 0.5, and a time increment
of 0.0005 seconds.

To compare the results in Reference 13 with the
present ones, the tip displacement of the beam in the Z
direction is plotted versus time in Figure 4—2. Either
graph in this figure corresponds to a value of B = 0.0625.
It is clearly seen that a good agreement exists between
the graph in the above—mentioned reference and the one
furnished by the present formulation (geometrically linear).

It is now of interest to study the axial force
effects on the response of the beam in both geometrically

linear and nonlinear cases. To this end, the

49

maximum absolute value of the tip displacements in the X

and Z directions are plotted against the axial load (static,

compressive) in Figure 4-3. Both graphs furnished by the

geometrically linear formulation correspond to a value of

8= 0. These two graphs are nearly straight while the a
I

graphs provided by the geometrically nonlinear formulation

are highly nonlinear, particularly, where the axial load

If the axial load exceeds 200 kips

 

approaches 200 kips.
(say, by 50 kips or more) the numerical results correspond-

ing to the latter formulation indicate that the beam would

collapse. It is to be noted that the axial load of 200

kips is equal to 40% of the axial carrying capacity (504
kips) of the beam cross-section and 26% of the Euler load

(771 kips, for the cantileVer case, of course) of the beam.

The numerical results based on the geometrically

linear formulation (not presented here) indicated that
the beam would collapse with an axial load between 400 and

500 kips. These loads are considerably larger than the

200 kip load indicated by the geometrically nonlinear

formulation. Thus, it may be concluded that when members

:arry substantial axial loads, the geometrically nonlinear

’ormulation, which is certainly more accurate, should be

sed.
To study the time-displacement response of the
the tip displacement in the axial (X) direction is

earn,
Both graphs plotted

-otted versus time in Figure 4-4.

50

correspond to an axial load of -30 kips as shown in Figure

In addition, the graph provided by the geometrically

4'10
This

linear formiation corresponds to a value of B = 0.

graph exhibits a permanent set of -0.002 inches which is

It also shows that there is no A

entirely due to yielding.
This

apparent physical vibration in the axial direction.

can be explained by noting that the axial shortening

 

effects do not enter the geometrically linear formulation.
(The high-frequency oscillations appearing in the above- ‘
mentioned graph are probably due to imperfections of the ._

numerical integration technique employed in the solution

procedure.) The graph furnished by the geometrically

nonlinear formulation shows a much larger permanent set

of -0.010 inches. This is partly due to the force inter-

action effects (-0.002 inches) and partly due to the axial

shortening effects considered in the geometrically non-

linear formulation.
The period of axial vibration predicted by the

geometrically nonlinear formulation is roughly 0.078

This period can be measured directly from the

seconds.
It is of some

relevant graph presented in Figure 4-4.
interest to compare this period with the one approximated

. 0
3y vaifiBIJEAx where ZmB is the mass of the beam

recording to the data chosen for the beam, the latter

>eriod is calculated to be 0.087 seconds. This period

lgrees fairly well with the one mentioned earlier.

51

inns comparison would thus mean that the beam appears to
respond.in the axial direction like a single-degree-of—
freedom system if the geometrically nonlinear formulation
is used.

It is recalled that in obtaining geometrically
nonlinear solutions, the axial force in a given time-
step, is approximated by the one obtained at the end of
the previous time-step. The question may be raised that
how much this approximation affects the accuracy of the
resulting solutions.

To consider such effects, a geometrically non-
linear solution of the beam was obtained without using
the above-mentioned approximation. This was accomplished
by performing iteration on the axial force within each
time-step. An axial load of 200 kips in compression was
applied as a part of the static loading in order to create
a more critical situation. The numerical results so ob-
tained (not shown here) did not indicate any noticeable
difference from those obtained with the approximation.
For example, the resulting displacements differed by less
Hunt 3% in the axial direction and 2% in the transverse
iirections. It is thus concluded that the performance
>f iteration on the axial force would not significantly
iffect the numerical results. The advantage of using the

Lpproximation, of course, lies in a significant saving of

:omputer time .

 

!‘ .
‘ I
' h

S Fifi
a

\

 

» . IA. ..
n _ .
‘qqguwnu.

z

 

52

4.2 Triangular Frame

 

The frame to be considered as a second example is
the one shown in Figure 4—5. The horizontal members
(girders) and the vertical members (columns) are fabricated
from 12WF53 and 12WF40 sections, respectively. The webs
of the girders are vertical. The webs of the columns lie Ei
in the vertical planes containing the bisectors of the i a

triangular planform. The external loads, in addition to

 

the weight of the girders, are shown in the figure. A
step-function type of dynamic loading is applied. The
value of a employed is the same as the one used in the
preceding example, that is, 0.5.

For the purpose of comparing the results in Refer-
ence 13 with the present ones, the displacement in the Z
direction of the joint 1 is plotted versus time in Figures
4-6 and 4-7. These two figures correSpond to the elliptic
and parabolic yield functions, respectively. The value
of 8 used in both cases is 0.125.

From Figure 4-6, it is seen that for the elliptic
yield function, a reasonably good agreement exists between
the two graphs. For the parabolic yield function, the two
graphs shown in Figure 4—7 are also quite similar. As far
as the displacement magnitudes are concerned, they differ
by approximately 7%.

In order to compare the geometrically linear and

nonlinear responses, the displacement of the joint 1 in

 

 

53

the Z direction is plotted against time in Figure 4—8.

The graphs presented in this figure correSpond to the
elliptic yield function. In addition, the graph furnished
by the geometrically linear formulation corresponds to a
value of B = 0 and a time increment of 0.001 seconds.

For the time increment just mentioned, the geometrically
nonlinear response turned out to be unstable. The largest
time increment for which the geometrically nonlinear
response became stable was 0.00005 seconds, i.e., 1/20
times 0.001 seconds.

It is seen that, the geometrically nonlinear
solution indicated a larger response than that given by
the geometrically linear solution. It should be remem—
bered, however, that the axial loads in this case are
rather moderate. As seen from the preceding example on
the cantilever beam, if the axial loads were sufficiently

large, the difference between the two solutions would be

much more drastic.

4.3 Building Frame

As a final example, a two-story two-bay building
frame, as shown in Figure 4-9, is considered. The member
9 (girder) of the frame has a 12WF53 cross-section. The
other girders and the columns are made from 8WF40 and 8WF17

sections, respectively. The webs of the columns are

marallel.to the south-north direction while the webs of

:he girders are vertical.

 

 

‘
'0 .1- TL
VF.— 1“

 

54

 

The static loading, in addition to the weight of
the girders, consists of a load of 1.50 kips per foot on
Umenember 9 and a load of 0.75 kips per foot on the rest
of the girders. The mass of these loads is lumped at
the relevant joints in the same manner as the mass of p
the members themselves. The dynamic loading is supplied he.
by subjecting the foundation of the frame to all three

components of the 1940 El Centro earthquake. The sketches

 

to be given in this example correspond to the elliptic .
yield function, a time increment of 0.0008 seconds, and
a value of a = 0.5.

Figure 4-10 shows a plan view of the distorted
shape of the frame as predicted by the geometrically
linear formulation. The solution corresponds to a value
of B = 0. In Figure 4-11 are shown similar results ob-
tained from a solution with the geometrically nonlinear
formulation. The distortions in both figures are based
on the member plastic displacements recorded at the end
of two seconds of the ground motion. It is seen that the
two distortions, suffered by the frame and computed by

the two formulations, are comparable as far as the general

appearance is concerned. However, as expected, the dis-

tortion magnitudes given by the geometrically nonlinear

formulation are greater than those of the linear version.

 

55

4.4 Cbmputation Time
In concluding the present chapter, it is instructive

‘U3cmnsider the computation time for the method given in

Reference 13 and that develOped in the present work. To

this end, for the three examples presented, the time incre—

ment, execution time, and real time interval of the several

solutions are listed in Table 4-1. (In all cases, the com-

puter execution time refers to the CDC 3600 digital com-

puter at Michigan State University.) It is seen from this

table that the execution time for the method presented in

Reference 13 depends on the value of 8, whereas for the

present method (geometrically linear), they do not. It

is also seen that the present method requires considerably

less computation time. This is perhaps the most signifi—

cant improvement of the present work over that reported

in Reference 13.

 

 

 

rifl

 

 

CHAPTER V
SUMMARY AND CONCLUSIONS

In this study, a matrix formulation has been pre—

sented for the dynamic analysis of space frames. In the

analysis, the effects of both geometric and material non-
linearities have been taken into account. A computer
program has been prepared for the implementation of the

analysis. Numerical results of three problems were ob—

tained in order to demonstrate the validity and practicality

of the formulation. These problems were: a cantilever

beam subjected to a pulse type of loading; a six-member
space frame subjected to a step-function type of loading;

and a two-bay two-story building frame subjected to the

1940 El Centro earthquake.

Comparative data with and without the effects of

geometric nonlinearities taken into consideration were

shown in.the form of graphs. Based on these graphs the

following observations were found to be noteworthy:
(1) The geometrically nonlinear formulation developed

in this work is obviously applicable to geometri—

cally linear problems. Good agreements were found

to exist between the numerical results obtained in

56

(2)

(3)

57

such applications and those reported in Reference

13 in which a different method was used. This may

be construed as an evidence for the validity of

the present analysis .

Plastic displacements as predicted by the geometri-
cally nonlinear formulation were generally larger

than those resulting from the geometrically linear

version. But, when axial loads were small, the

influence of geometric nonlinearities on the

numerical results presented herein did not seem

significant. However, as axial loads increased,

the influence rapidly grew. For the cantilever

beam problem, the beam would collapse with an

axial load equal to approximately 26% of the Euler

load. This axial load is practically equal to

only one half of the magnitude corresponding to

a geometrically linear solution.

The numerical results of the cantilever beam problem

based on using the axial force of the previous time-

step, for the calculation of the member stiffness

matrices, differed only insignificantly from those

resulting from the more accurate approach of iterat-

ing on the axial force. The advantage of not per-

forming any iteration lies, of course, in a con-

siderable saving of computer time.

 

L

 

 

58

Natural extensions of this work may include the
incoquation of more accurate yield surface equations in
thecxmputer program, and an investigation into the possi—
kfiJJty of reducing the number of degrees of freedom in
order to facilitate applications to even larger structural

systems such as high-rise building frames.

 

 

 

59

Tabhe4-l Comparison of Computer Time between Reference 13
and Present Work

 

Computer . "Real Time"
. . Time Increment .
Execution Time . of Solution
. in Seconds .
in Seconds in Seconds

 

Cantilever Beam

 

 

Reference 13 (B=l/8) 11 0.0004 0.06
Reference 13 (B=3/32) 14 0.0004 0.06
Reference 13 (B=l/l6) 170 0.0004 0.06
Present Work 4 0.0005 0.08

Triangular Frame
Reference 13 (8:1/8) 153 0.0005 0.4

Present Work 23 0.001 0.

Building Frame
Reference 13 (8:1/8) 998 0.0005 2.0

>resent Work 315 0.0008 2.0

 

 

60

 

 

 

:{> X

Joint Coordinate Axes

Z 64) 3 Member M

J : Joint J

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Member Coordinate Axes for Member 3

Figure 1-1 Typical Space Frame and Notation

61

 

 

 

 

 

. y2
Y2 A
\ x'
\ ’57 2
\ //
/
yl £>x2
A / °2
/
/
22 1!
22 Y
#D x
01 1
z
/ "

"51

Figure 1-2 Different Coordinate Systems

 

4
ball!

 

62

- . ...A. ...p
.... Iltma.E .
F

O
H

I

.

uzs

 

 

 

x2:

 

 

NW9

 

 

 

 

 

um3.Nm.z\Q

 

 

 

 

63

 

 

 

 

 

 

 

 

 

Y
A
(l+e)dx
A B
——(> dx <——— > x
dx <——— u M D X
I A B x V4

 

 

 

 

 

u<:

 

 

Figure 2-2 Free-Body of a Differential Element

 

 

 

 

 

 

64

 

 

M
A y'”)!
Mx = constant
M = constant
2 AF / V
AUP
F
(>
F , u
x x
J

 

Yield Surface 4 a 1

Figure 2-3 Geometrical Interpretation of Flow Law

65

me mocmnomcm cw pomp mm Hopoz oumuomwo vim ousmwm

 

 

 

 

 

 

    

mumm cease muaaanaxoam msosaapaoo

cues coHuwom owummam

 

suaannaxoam ceases sues coauuom oaummam

66

 

b ucHOH cues pouoHOOmmm moon cwmflm Hum ounces

Houcou who: b boson

 

 

 

 

 

-;Hil
meadow mum: m cocmum c
n I
.B AB
R n N
AB N
in» the ..noc n he s xA
new Ln». Luvs n he

 

w«<»

67

Yield Surface 4): l

l

\

2)

 

 

igure 3-2 Treatment of Force Vector during Yielding

 

 

68

 

 

 

2 12WF53 [E]

 

 

ax
A(//////VGlobal
Coordinate

Z System

UJ.__L.Lb
V

 

em- .. eve—4

 

1. .-
i

I
Joint Axis Static

 

 

Dynamic Duration of
Load Load Dynamic Load
1 X -30 kips -- --
1 Y 1 kip -18 kips 0 to 0.01 sec.
1 Z l kip 8 kips - 0 to 0.01 sec.

 

 

 

 

 

Figure 4-1 Cantilever Beam and Loading

69

 

mafia mcmuo> Econ Ho>oawucmo mo ucchochmflo omno>mccus mlv madman

mpcoomm as mafia

mo.o v0.0 No.o

 

 

_ 14 fl

scones haacofluuosomw

.m whamflm Scum (llIIqu

mmmo.o u a .ma oncogenes

\
I‘ll.
tl‘ lil‘

 

 

«tn

seqou]: ujt uoraoerra z u; queuteoetdsra dm,

70

 

 

 

  

    

     

 

 

 

 

I I — II.
r. r
a a
e e
n n
.1. .1.
L L
m m L
a V. a
C 1 C
.i .1 .i
Y r a r
n e .... u ....
a m r
.m a e um m ..
r e G 1
+.n m." c
e .1. e O
mflm a.N
e O
G N /
8 4 A0 nw nw AN 0
0 0 a. 9. l
O O
0 GA HGOEGOMHQ
O coauomwdo N .
m0 OQH d.“ .
monocH . . . n4 Eocexmz .

(men mfle

m0 05Hm> QUSHOm

100 150 200

50
Axial Load (Static, Compressive) in Kips

Figure 4-3 Effects of Axial Load

 

mafia msmuo> Econ um>oawucmo mo ucosoooammwo Hmflxc «Iv ounces

 

71

 

 

 

moonwacoz asamoauumeooo .IIIIq

>15 A}

imoo.ol

 

wooden
mascoauuchou.llJ

\ ..a R .

 

 

lilitLu‘ . o
.Zlivhllllilllllllllllii V O o O NO 0

mpcouom cw mafia

smpuI u; uorioerra x u; quemeoetdsrq d1;

 

72

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Y
$1 Global 2
Coordinate
System
j> x
g/ ) 60 60
Z / I 3
5
11ijan 10,
'<cs 14' a:
FWHHTn nTHTn1
4 6
Joint Axis Static Dynamic Duration of
Load Load Dynamic Load
1 Y -30 kips -- --
1 z 1 kip 20 kips 0 to w sec.
1 Y 6 kip-in. 120 kip-in. 0 to a sec.
2 Y -30 kips -- --
3 Y -30 kips -- --

 

 

 

 

 

 

Figure 4-5

Triangular Frame and Loading

 

 

 

 

4.» rflrrlbwflm N."

 

Joint 1 Displacement in Z Direction in Inches

2.8

2.4

73

 

—

 

 

B = 0.125

Geometrically Linear

 

'-"—"' From Figure 11, Reference 13

 

 

l I l
0.2 0.3 0.4

Time in seconds

Figure 4-6 Response of Triangular Frame

(Elliptic Yield Function)

 

 

(

 

Joint 1 Displacement in z Direction in Inches

 

 

2.8

2.4

2.0

1.6

1.2

 

 

B = 0.125
Geometrically Linear

From Figure 11, Reference 13

I l l l
0.1 0.2 0.3 0.4

Time in seconds

Figure 4-7 Response of Triangular Frame

(Parabolic Yield Function)

 

 

JOint 1 Displacement in Z Direction in Inches

75

 

 

 

 

 

Geometrically Nonlinear
2.8t' //
2.4“-

200 —
1.6 _
Geometrically Linear
(8 - 0)
1.2 —
0.8 —
0.4
0 J l 1 l
0.1 0.2 0.3 0.4

Time in seconds

Figure 4-8 Comparison betwaen Geometrically Linear

and Nonlinear Responses of Triangular Frame

 

 

 

76

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Y’ UP / North
/'
/
//
West <——_| —> X, East
| Global
$ Coordinate
System ..
Down
Z, 80111111 2 3
1
1
I
1 12
1. ., g
o 6
14'
9 5
.4". men .J
20'
”(m .1 A/

 

 

 

I<e 20'

Figure 4-9 Building Frame

 

h'n nT'rn
20'—fl

 

 

77

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(4.27,1.16) (3.11,0.70)
r—DX :- ————————————— '1
I 2 F i
2 I 3
1 Plastic ' 2
' Displacements I 1
Shown in . I
Parentheses I I
z in Inches !
i I
I l
1 i i 4
L.___._.JL ________ _.l
4.60 0.23
( ' ’ (3.35,o.32)
,_(§_-_OZ_r_°.-_5_51 ______ TC} 2.7143423). _____ 121- 94 , 0 . 56)
.’ I I
I
. E: :7
I l I
l I I
I l I
I l l
I I I
l I l
I I I
I I I
I I l
E]: 5 I i 8
L_.__ a..- ________ .l.__.__ _________ _J
(2.92,0.62) (3.97,0.72) (3.35,0.32)

Figure 4-10 Plan View of Distorted Frame Predicted
by Geometrically Linear Formulation

 

 

 

I... ______________

 

 

 

 

 

2

 

 

3

 

Permanent

 

Displacements
Shown in

 

Parentheses

78

x (5049,1005) (5005,0072)
I
l
|
l
I
l

in Inches I r
l I
I
I
l
I
I

 

 

—-———.~—u— —-——————. F.
- 1..-. t ‘n' 'T
C- a
._ ...
- a

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Z
i
'1 4 ..
L_____i ________ J -i
(5.69,0.53) (6.38,0.37) 3
(2:393:19) ________ (52.21295). ______ <5-06r0-68I
I" I ‘I
l 5 I l 7
I I l
I : I
I I
I I l
I
I ' I
I ' I
I ' l
I ' I
I 5 I I
L.-- _ __________ .L _______________ J 8
(2.82,0.55) (5.28,0.95) (6.38,0.37)

Figure 4-11 Plan View of Distorted Frame Predicted
by Geometrically Nonlinear Formulation

 

 

LIST OF REFERENCES

 

f3 .

-/ 4 .

LIST OF REFERENCES

Argyris, J. H., Recent Advances in Matrix Methods of

Structural Analysis, Pergamon Press, The Macmillan
Company, New York, 1964, pp. 115-145.

Connor, J. J., Logcher, R. D., and Chan, S. C.,
"Nonlinear Analysis of Elastic Framed Structures,"
Journal of the Structural Division, ASCE, Vol.

94, No. 8T6, Proc. Paper 6011, June, 1968, pp.
1525-1547.

Hodge, P. G., Plastic Analysis of Structures, McGraw-
Hill Book Co., Inc., New York, 1969.

Iverson, J. K., "Dynamic Analysis of Nonlinear Elastic

Frames," Ph.D. Thesis, Department of Civil Engineer-
ing, Michigan State University, 1968.

Jennings, A.,"Frame Analysis Including Change of
Geometry," Journal of the Structural Division,

ASCE, Vol. 94, No. 8T8, Proc. Paper 5839, March,
1968, pp. 627-644.

Johnson, D., and Brotton,[L NL, "A Finite Deflection
Analysis for Space Structures," International
Conference ong§pace Structures, Department of
Civil Engineering, University of Surry, Surry,
England, September, 1966.

 

Livesley, R. K., Matrix Methods of Structural Analysis,

Pergamon Press, The Macmillan Company, New York,
1964, pp. 230-258.

Malvern, L. E., Introduction to the Mechanics of a
Continuous Medium, Prentice-Hall, Inc., Englewood
Cliffs, New Jersey, 1969, pp. 356-363.

Morris, G. A., and Fenves, S. J., "A General Procedure
for the Analysis of Elastic and Plastic Frameworks,"
Structural Research Series No. 325, Department of

CIVil Engineering, University of Illinois, August,
1967.

79

WT.
!

 

 

10.

11.

12.

‘I 13.

14.

15.

80

Morris, G. A., and Fenves, S. J., "Approximate Yield
Surface Equations," Journal of the Engineering
Mechanics Division, ASCE, Vol. 95, No. EMZT Proc.
Paper 6741, August, 1969, pp. 937-954.

Nigam, N. C., "Yielding in Framed Structures under
Dynamic Loads," Meeting Preprint 737, ASCE Annual
Meeting and National Meeting on Structura1 Engineer-
ing, Pittsburgh, Pa., September 30-October 4, 1968.

Saafan, S. A., "Nonlinear Behavior of Structural Plane
Frames," Journal of the Structural Division, ASCE,
Vol. 89, No. 8T4, Proc. Paper 3615, August, 1963,
pp. 557-580.

Wen, R. K., "Dynamic Response of Elasto—Inelastic
Structures," Progress Report No. 1, NSF GK-590,
Division of Engineering Research, Michigan State
University, February, 1969.

Wen, R. K., and Iverson, J. K., "Space Frame - Statics,"
User's Manual, Department of Civil Engineering,
Michigan State University, February, 1969.

Zarghamee, M. S., and Shah, J. M., "Stability of Space-
frames," Journal of the Engineering Mechanics
Division, ASCE, Vol. 94, No. EMZ, Proc. Paper 5878,
April, 1968, pp. 371—384.

 

 

‘5. w. s . A211Jefii-
‘
‘1‘ 7112-“.

APPENDICES

 

j".

 

I
1 ~ I _
‘ L}‘I_t-'J' ‘ I

 

APPENDIX I
STIFFNESS MATRICES

In this appendix, the elastic stiffness matrices

are expanded, and the geometric stiffness matrix is de-

rived.

A.1.1 Elastic Stiffness Matrices

The terms employed in the elastic stiffness matrices

S etc., are first defined. The matrices themselves are

PP’
expanded next. The expanded form of the matrix FG is also

given.
3 = Y: Z
k = y, z (A-l—l)
3'7! k
i. 2" 0.5
_ I FNXL 1 2
(1)3.- L— T (A- ")
C. == c s .
3 0 (DJ
(A—1—3)
S. == sin ¢-
3 J

81

 

 

‘“ :‘l.-
-‘\

..__

w 1

his-I":

Er.

 

82

¢-(¢. - S ) .
hjy '3 fi-ijag-yT-g-T (3-1-4)

:3"
II
{T
+
5'

(A—l-S)

(py and pz are the member chord rotations in the x—z and

x-y planes, respectively.)

(A—l-S)

 

 

 

 

 

 

 

 

t
FG [EAx(6y+6z), oz FNX, py FPX, o, o, 0] (A-l 7)
EAx;
L ... L. ...——
‘ Eiz EI
2h 2 -_5 h -§
Y’ L ”fly” _ yz L
EI EI
S :- thz L -1122 L
PP
cg
L
EI
11 .;JL
zx L
symmetric h BI;
yx L

 

 

 

 

 

 

(A-l-B)

 

s~'__' [IFI_

 

 

83

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

-EA3
"E"
312 EIz
“In? “an?
El E1
”hm—55 ”22—55
SPN ' -— _GJ
..---_ f-
I E1 E1
I hzz L! hzy L a
..h EIz EIz :
2“?” L
L. Y L I yy _J
(A-l-9)
L ; I
1 ”Wt—7-7 .1 .1 .4 .- ...--- -
312 I 312
“w? ' ‘52?
I E1 1‘ ’EI ””””“‘”
I
i2h 1 h -—%§
V i 22 L 1 22 L I
SNN " ' elf ‘
1 1-._1-1- 11-11-...L11I1_-_._-...--
I EI
---..- I hzx L
I EIz
symmetric I hYXT
_, I ._
(A-l—lO)

IA.1.2 Geometric Stiffness Matrix

To arrive at the geometric stiffness matrix sG
used in this study, the incremental rotations due to

flexure are neglected, that is,

84

Am = A
y py

(A-l-ll)

sz = Apz

The incremental chord rotations Apy and Apz are constant
along their corresponding member. The approximation men-
tioned above is consistent with neglecting the change in k.
the elastic stiffness matrices SPP’ etc., when forming the

increments.

 

Consider now the matrix PG produced in the previous

section. Its increment is given by

AFG = [EAX (A6y + A62).

Px Px'
t
0, 0, O] (A-1-12)
in which
_ l L
Aéy - '11: 6 (Dz sz dX
Apz IL du
= _L— 0 dx
= 0z ADz
O
z (A—1—13a)

u
I
r'I
ST
5
I<

I

D
2?

<V

P

and similarly,

A6 =
z

and also,

AF =

PX

and

AFNX =

It is now easy

into

85

D
‘LXI(AuPz - AUNz)

EAX
L (Aqu - AuNx)

- EA (A6 + A6 )
x y z

EAx
L (Aqu _ AuNx)
EAx
+ L (AuP Au‘Ny) pz
EAx
- L (mfiz-A%M)py
-AF

PX

to show that the increment AF

-S (AUP - AUN)

(A-l-l3b)

(A-1-14a)

(A-l-14b)

can be cast

(A-1-15)

(A-1-16)

'1 -' “Lorr- Inn-1'
‘ .

- 71"

 

qr-

.—..__._

where

86

To

symmetric

 

OZ

2 + FNx

02 EA
X

 

(A-1—17)

. IV! ‘.1‘h.'_bl=‘n

 

APPENDIX II

INCREMENTAL FORCE-DISPLACEMENT RELATIONS
FOR A GEOMETRICALLY MEMBER

WITH END RELEASES

‘2‘

Assume that a geometrically nonlinear member PN,

‘3“:
‘.

as shown in Figure 2-1, has a total of R releases at its
negative end. Corresponding to the Ith release, assume

a column matrix V defined by

I

vI = [0, o, o, 1, o, 0]t (A-Z-l)

in which the only nonzero element is a "1" corresponding
to the released force component (at the negative end of
the member). For instance, the matrix VI shown above
corresponds to a released torsional moment.

The negative-end displacement increment AUN can
be decomposed into an elastic part Au; and a part AU;

due to the releases being considered. Thus

e _ _ r - -
AUN — AUN AUN (A 2 2)

The increment AU; may be written as

R
= 2 v b (A-2-3)

87

88

in which bI is a scalar factor which gives the magnitude
of the displacement increment at the Ith release. It is

now convenient to define a 6xR matrix G by
,v ] (A—2-4)

If the factors bI (I=l, ..., R) are assembled into a k~

column matrix B, Equation (A—2—3) reduces to

AU; = GB (A-2-5)

 

Substitution of Equation (A-2-5) into Equation (A-2-2)

furnishes

e- - ——
AUN — AUN GB (A 2 6)

The latter equation is now introduced into the incremental
force-displacement relations given by Equations (2-17).

This leads to

AFP = KPP AUP + KPN (AUN-GB) (A—2-7a)
- t _ _ -
AFN - KPN AUP + KNN (AUN GB) (A 2 7b)

The condition that the negative-end force component
increment in the direction of any release I is zero is

expressed by

VI AFN = 0 (A-2-8)

89

Substituting Equation (A-2-7b) into Equation (A-2—8) and

considering all the releases, it follows that

 

t t t _
VI KPN AUP + VI KNN (AUN-GB) — O
I = 1, ..., R (A-2-9)
These equations can be combined into the following single
equation:
t t t _ - - -
G KPN AUP + G KNN (AUN GB) — 0 (A 2 10) E;
from which
_ t t t - -
B — E G KPN AUP + E G KNN AUN (A 2 11)
where
_ t -l - _ _
E - (G KNN G) (A 2 12)

Finally, substituting Equation (A—Z-ll) back into Equations
(A—2-7) , the following incremental force-displacement re-

lati ons are obtained:

AFP = (KPP - KPN G E G KPN) AUP
+ (K - K G E Gt K ) AU
PN PN NN N
(A-2-l3)
_ t _ t t
AFBq — (KPN KNN G E G KPN) AUP
+ (K - K G E Gt K ) AU
NN NN NN N

a
”any“

 

 

II..__..sI_-» . -_-I

 

 

 

 

90

The extension to the case where the positive end
of the member has also releases is straightforward and

thus omitted.

‘_.— 2:. m @

s-ufi'r.

 

“:31!!!-

APPENDIX III

COMPUTER PROGRAM

For the sake of completeness, the computer program
written for the present study is outlined in this appendix.

'Dhe routines which constitute the program are briefly

described. The important identifiers used in the program

are defined; and finally, a listing of the program is
presented.

A.3.1 Description of Routines

The computer program developed consists of a main

routine called DYNAMIC and four subroutines named IN-

VERSE, RMATRIX, FMATRIX, and SMATRIX. The first sub-
routine INVERSE, taken from the M.S.U. computer laboratory,
is used to invert square matrices as needed in the program.
It is the only subroutine written in COMPASS (CDC 3600
machine language); therefore, its listing will be omitted.

The second subroutine RMATRIX is used to compute the

responsiveness matrix for a given joint.

The third subroutine FMATRIX has two versions.
One of the two corresponds to the geometrically linear
fonmfletion (of a discrete model as described in Section

2.3)vmile the other is concerned with the geometrically

91

 

92

nonlinear formulation. The linear version is used to

compute the stiffness matrix for a given member. Note

that the eight nonzero elements of the latter matrix are

the only ones which are actually computed. The non—

linear version of the subroutine is used to compute cer-

tain factors (a total of eight) that appears in the elements L

of the elastic stiffness matrices for a given member.

These factors are independent of the axial force. The

effects of the latter force are considered in the sub-

 

routine SMATRIX.

Analogous to the subroutine FMATRIX, the fourth

subroutine SMATRIX has also two versions. Both linear and

nonlinear versions of this subroutine are used to perform
the same task, namely, to compute the member end forces.

They are also capable of inserting and/or removing plastic

hinges when necessary. Each version of this subroutine

has in turn two versions depending on the (type of) yield

function to be used. In the listing, to be presented

shortly, the linear version of the subroutine SMATRIX
corresponds to the parabolic yield function while the

nonlinear version corresponds to the elliptic yield func-

tion.

In the main routine DYNAMIC, the initial static

sohnfion is first read in as a part of the input. This

sohndon is conveniently obtained by using the program

STATIC written by Wen and Iverson (l4) . The rest of the

93

input is read in next and immediately printed out. The

identifiers are then initialized as necessary. In the

following portion, the output of the program is monitored

for printing and plotting. If the integration time exceeds

the time limit specified, the program stOps. The incre-

mental joint displacements are computed next. The rest L

«of the program.is concerned with computing: the joint

forces due to the member end forces, the ground accelera-

 

tion if applicable, and the joint accelerations and veloci-

ties. The control is then sent back to the portion where

the output is monitored.

A.3.2 Definition of
Important Identifiers

The important identifiers used in each routine are

defined below in alphabetical order. Any identifier which

appears in more than one routine but maintains the same

meaning is defined only once.

PROGRAM DYNAMIC

ACC(I,J) = Ith component of joint J acceleration

ALPHA = a

AMPLIFY = amplification factor times acceleration
of gravity

AREAX(M) = cross-sectional area of member M

AREAYUH =

effective shear area of member M in y

direction

AREAZ (M)

BETA

COMP (I ,M)

DATA A(I,J)

DATA T(I,J)
DDIS(I,J)

DEADM (M)

DENSITY
DIS(I,J)
DPLT
DPRT

DT

E

FNN (I ,M)

FORCE(I,J)

FPP(I,M)

G

GLINEAR

94

effective shear area of member M in

z direction

8

projection of member M on Ith global
axis

Ith component of foundation acceleration
in Jth reading

time correSponding to DATA A(I,J)
increment of DIS(I,J)

uniformly distributed mass per unit
length of member M

density of material

Ith component of joint J displacement
time interval for plotting

time interval for printing

time increment

modulus of elasticity

Ith component of force on negative end
of member M

Ith component of force on joint J
exerted by incident members

Ith component of force on positive end
of member M

modulus of shear

true, if geometrically linear response

desired; false, otherwise

 

95

IXX(M) = torsional constant for member M
cross-section

IYY(M) = moment of inertia about y axis for
member M cross-section

IZZ(M) = moment of inertia about 2 axis for

 

member M cross-section

 

JFREE = number of free joints

JN(M) = number, identifying negative joint of
member M

JOINTS = number of joints

JP(M) = number, identifying positive joint of
member M

LENGTH(M) = length of member M

MEMBER(I,J) = number, identifying Ith incident member

on joint J; it is signed according to

incidence
MEMBERS = number of members
MLINEAR = true, if material nonlinearity effects

ignored; false, otherwise

NMEM(J) = number of memebers incident on joint J

PEXT(I,J) = Ith component of dynamic load on joint J

QUAKE = true, if earthquake applied; false,
otherwise

QTIME = time lag

RESPON(I,J,K) = element in Ith row and Jth column of

joint K responsiveness matrix

 

 

RM(I,J,M)

ROT(1)

ROT(2)

SLIMIT(I,M)

STLOAD(I,J)
T(K)

TIME
TLIMIT
VEL(I,J)
xo ACC(I)
x0 DIS(I)

X0 VEL(I)

96

element in Ith row and Jth column of
rotation matrix (3x3) associated with
member M

cose, in which 6 = angle measured clock-
wise from east—west axis to global X axis
sine (see above for definition of 6)

Ith fully plastic stress resultant for
member M cross-section

Ith component of static load on joint J
Kth component of foundation acceleration
time

time limit

Ith component of joint J velocity

Ith component of foundation acceleration
Ith component of foundation displacement

Ith component of foundation velocity

SUBROUTINE SMATRIX

 

(linear version)

33(1)

33(2)

CHI(1)

plastic-flow constant corresponding to:
either section x = BL, or section x =
(l —B)L if section x = BL not yielding
plastic-flow constant corresponding to
section x = (1 - B)L if section x = BL
yielding also

yield—function value for force on

section x = 8L

 

“ir

 

WK

 

CHI (2)

ITV(I)

FP (I)

GG(I ,J)

ITE RATE

LOOP

PDISP(I,M)

PHI (l,M)

PHI(2,M)

PLASTIC(1,M)

PLASTIC(2,M)

sum

SS (I,J)

V(I,l)

97

yield-function value for force on section
x = (1 - B)L

Ith component of member distortion

Ith component of positive-end force
element in Ith row and Jth column of
matrix defined by Equation (2-37)
false, initially; true, if member end
force not acceptable

identifier to prevent endless lOOping
Ith component of plastic distortion of
member M

CHI(1), if section x = BL of member M
not yielding; constant, otherwise
CHI(2), if section x = (l - 8)L of
member M not yielding; constant, other—
wise

true, if section x = BL of member M
yielding; false, otherwise

true, if section x = (l - B)L of
member M yielding; false, otherwise

Ith stiffness element for member M
element in Ith row and Jth column of
member stiffness matrix

Ith component of yield surface normal
corresponding to: either section x = BL,
or section x = (l - B)L if section x =

BL not yielding

 

[Pfi_ nu“
(

“we.

 

h”?

 

V(I,2)

YIELDED(1,M)

SUBROUTINE SMATRIX

98

Ith component of yield surface normal
corresponding to section x = (1 - B)L
if section x = BL yielding also

false, before first yielding occurs in

 

member M; true, afterwards L

 

(nonlinear version) p

 

BB(J)

CHI(J)

DUN(I)

DUP(I)

ENDN

ENDP

ENDS

FN(I)

GP(I,J)

plastic-flow constant; J = 1 corresponds g

 

to positive-end cross-section, J = 2 34‘
corresponds to negative-end cross-section
yield-function value; J = l correSponds

to positive-end force, J = 2 corresponds

to negative-end force

Ith component of negative—end diSplace—

ment increment

Iflicomponent of positive-end displacement

 

increment

true, if only negative-end section yield-
ing; false, otherwise

true, if only positive—end section yield—
ing; false, otherwise

true, if both end sections yielding;
false, otherwise

Ith component of negative-end force
element in Ith row and Jth column of

matrix defined by Equations (2-26)

| 11' '1' l'. Eli. .I.) II "Illl lull! l I. III

GN(I,J)

H21, ..., H33

P

PDISN(I,M)

PDISP(I,M)

RHOZ, IdHMB

S(I,M)

SNN(I,J)

SPN(I,J)

SPP(I,J)

VN(I)

VP(I)

YIELDED(l,M)

99

element in Ith row and Jth column of
matrix defined by Equations (2-26)
stiffness factors

axial force

Ith component of negative-end plastic
displacement of member M

Ith component of positive-end plastic
displacement of member M

py and pz, respectively

Ith common factor of elements of elastic
stiffness matrices for member M
element in Ith row and Jth column of
direct stiffness matrix, referred to
negative end

element in Ith row and Jth column of
cross stiffness matrix,referred to
positive end

element in Ith row and Jth column of
direct stiffness matrix, referred to
positive end

Ith componenent of yield surface normal
corresponding to negative—end section
Ith component of yield surface normal
corresponding to positive-end section
false, before first yielding occurs in
positive-end section of member M; true,

afterwards

 

 

 

lOO

YIELmMM2,M) = false, before first yielding occurs in
negative-end section of member M; true,

afterwards

SUBROUTINE RMATRIX

 

BIM(K,L) = element in Kth row and Lth column of
branch inertia matrix

JIM(K,L) = element in Kth row and Lth column of
joint inertia matrix

MMASS(I) = mass of branch I

UX, UY, UZ = coordinates of joint mass center with
respect to local coordinate system
through branch mass center

WX, WY, WZ = coordinates of joint mass center with

respect to joint coordinate system

 

 

 

101

A. 3.3 Listing of Prggram

PROGRAM DYNAMIC
LOGICAL QUAKE. GLINEAQO -MLINEAR9 GNONLIN
DIMENSION ACC(6050)0 VEL(6950). PEXTIécSO)! FORCEI60‘30)

DIMENSION DATA T(39800)o DATA A(39800)o ROT(2)9 T(B)
DIMENSION XO ACC(3)9 X0 VEL(3). XO DIQIZBIQ INDEXIB)
DATA (X0 DIS = 3(0)). (XO VFL = 3(“))

DATA (X0 ACC = 3(0)). (INDEX = 3(1))

COMMON REQPONIéQéoQO)

COMMON/ TIME/ TIMF’Q DT

COMMON/ STIFF/ 9(Rolnfi)o B(?OIOO)
COMMON/ PLASTIC/ PLAGTICIZQIOO) $ LOGICAL PLASTIC

COMMON/ YIELDFD/ YIELDEDI29100) 5% LOGICAL YIELDFD

COMMON/ DIS/ DIS(fio‘SO). DDIQIGQSO).
I pDI§p(6QIOOIQ pDIgNIGOIOO’

C I
c IDENTIFIERS WHOSE VALUE: TO BE READ IN
c FROM INITIAL STATIC SOLUTION

0

COMMON/ TYPE/ DENSITY. E. G
COMMON/ QI7E/ MFMRFRQQ JOINTfio JFRFF
COMMON/ ARFA/ ARFAXI 100) 9 ARFAYI 100). ARFAZI 100) o

I Ixxunn). IYY(IOO)9 I77I100) $ QFAL IXX. IVY. I77

COMMON/ DFADM/ DF'ADMIIOO)
COMMON/ COMP/ COMDCBQIOD)
COMMON/ SLIMIT/ QoLIMI'IWFnIOO)
COMMON/ LENGTH/ LENGTH(IOO) 5 PEAL LENGTH
COMMON/ FORCE/ FPPIfioIOOIv FNNI6OIOO)
COMMON/ JPJN/ Jp(100)0 JNIIOO)
COMMON/ DM/ RM(303QIOO)
COMMON/ MFM/ NMEMCRO)Q MFMRFPIIOQRO)
COMMON/ STLOAD/ 9TLOADI6Q'3O.)
C—--- -------------- C
C PDOGQAM INPUT C
C ————————————————— C
LUN = '3‘;
CALL SKIPF’ILEI LUN )
CALL SKIPFILEI LUN I
READ(LUN) MEMBERS. JOINTQO JFREE. DENQITYQ E9 6
DO 505 M 2: IOMFMBEPS
505 READCLUN) AREAXUW). AREAYIM). AREAZHMI). IXX(M)9 IYY(M)0
IZZ(M). DFADMIM)9 (COMPIIQM). I=193)9
(SLIMIT(I0M)0 I=lo6)o LFNGTHUW).
(FPPIIQM)0 I=106)9 (FNNIIQM)9 I=I.6)9
JPN“). JN(M). ((RMIIQJQM)9 1:193). J=1g3)
DO :15 J = IQJOINTS
=51; PFAD¢LUNI NMEM(J)9 (MEMRERIIQJIQ I=I¢IOIQ
I (STLOADIIQJ)9 I=106)

 

 

“:._:‘."flu ~ '-".

_Jl

102

READ 531. QUAKE. OTIME.
READ 531. GLINFAR. BETA.
READ 531. MLINEAR. ALPHA.
FORMATILIO.BEIO.5.AOXI
IFI QUAKE I READILUNI
PRINT 537. QUAKE. OTIMF.
I CLINEAP. BETA.

2 , MLINEAP. ALPHA.
FORMAT(*IQUAKF=*L?*
I*OGLINEAR=*L2* HETA=*F7.4*

2*OMLINEAR=*L2*

511
POT.

537

D $ DT DTIME
PLTIME =-O.fi*DTIME
ONOTOGLINEAP
I.MFMBEP§
BETA*LENGTH(MI
LENGTH(MI BII.MI
PLASTICII.MI PLAQTICI?.MI
YIELDED(I.MI YIELDED(?.MI
CALL FMATRIXI M )
DO 640 I 1.6
PDI5P(I.MI PDIgNII.MI
SLIMITII.MI

PPTIME

GNONLIN
DO 640 M
BII.MI
8(2.MI

O

640
C
C
C

EXTERNAL DYNAMIC LOADS AT TIME

I.JFREE $ DO 670
0

DO 670 J
PEXT(I.J)

670
C
C
C

JOINT ACCOO VELOO AND DIQQ

DO 680 J
DIQII.JI
EN) 690 J
CALL PMATR I XI
00 690 I 1.6
ACCII..H VEL(I.JI
[X3 685 K 1.6

ACCII..H ACC(I.JI
IF( I 0LT. 4 I ACCII.JI
GCMI CCWTTINUF

C--wn- -------------- C

C DDOGDAM OUTPUT
C ------------------ C

I_Uhl PR

W9 I TE (LUN. 2005) JFREE
IFI 'TIME .LT. PRTIME I

n

I.JFPEE
J. ALPHA I

69O

O

6R=i

Ponn GO

AMPLIFY
DPRT.
DTIME.

AMPLIFY.
DPRT.

DTIME.
QTIME=*F704*

ALOHA=*Fq o 2*

AT TIME

I.JOINT§ Q DO 690 I

DPLT

TLIMIT

XO ACC. DATA T. DATA A

DPLT.

TLIMIT

AMPLIFY=*F7.2/

DPLT=*F8.5/
TLIMIT=*F8.5I

 

DPPT=*F8.5*
DTIMF=*F8.5*

 

OPALQF. o
OFALQF.

IoO/SLIMITCI.MI

I

+ RESPONII.K.JI*PEXT(K.JI

ACCII.J) - XO ACC(II

TO 203m

 

 

.: lull-

 

103

PQTIME - OPTIMFI + DPDT

RTHW'= TIMEF¢AI

PRHW'?OOI. TIME. RTIME
FOWW"W*IDYNAMIC VAQIAHLFS AT TIME =*F8.5.75X

2001

1 *ELAPSED TIME =*-3PF7.2I

IF( QUAKE ) PRINT 2002. X0 DIS. X0 VEL. XO ACC
20mg FORMAT(*OGROUND MOTION*30X*X*15X*Y*15X*Z*//

l I8X*DISPLACFMFNT*RX3F16.8/

2 18X*VELOCITY*1PX3F16.R/

3 19X*ACCFLERATION*RX3F16.R)

POINT 9003
fififiq FORMATI*-FORCF* 5x*<1*14x*<?*14x*<1*14x

*sa*14x*ss*14x*86*14X*Sl1*11X*SI?*)

(M. (FPP(I.MI. I=I.6I.
FNNI5.MI. FNN(6.MI. M=I.MFMRERSI

1
PR INT 2004 .

I
2004 FORMAT(/I4.BF16.8I
2008 FORMATI I4.6F16.8I

C
PRINT 2006

?006 FORMAT(*-DISPLACEMENT*I
PRINT 2007. (J. (DIS(I.JI. I=l.6). J=I.JFREFI

2007 FOPMAT(/I4.6F16.BI

C
IF( GLINEAR I PRINT ?012

IF( GNONLIN I PRINT 2013
FORMAT(*-PLASTIC DIQTORTION*I
FORMAT(*-PLASTIC DIQPLACFMENT*/I
DO 2016 M = I.MEMBER9

L =-M
IF( YIELDED(1.M) I PRINT 2004. M. (PDISP(I.M).
IF< YIELDED(2.M) I PRINT 2005. L. (PD!SN(I.M). I=l~6)

?016 CONTINUE

C
PO30 IF( TIME .LT. PLTIMF I GO TO 2040

20]?
ZOIB

1:1.6I

PLTIME = DLTIME + DPLT

WDITEKLUN.203PI TIME. ((DISII.JI. I=1.6I. J=19JFREEI
203? FORMAT(F'R.6. 6FI?.R/(8X 6FI?.8II
2040 IF( TIME .GE. TLIMIT I RETURN
c ---------------------------------------- r
C I NCREMENTAL FREE-JO I NT 0 ISPLACEMENTS C
C ---------------------------------- - ----- F

1001 J = I.JFREE $ DO 1001 I = 1.6
I*DT

DO
IFMDI IDDISS(I.JI = (
TIA": = TIME + DT

c—-—---—---- --------------------------------- c
c Foss—JOINT FORCES our- TO MEM. FND FORCES c

c—--_—-———-- --------------------------------- r
CALJ. SMATRIX( BETA . MLINEAR. GLINFAR I

Do 1506 J = 1.JF’REF

on 1501 I =1.6

VEL(I.J) + 0.s*ACC(I.J)*DT

 

104

I‘mI FORCE(I.JI = 0
Fl: NMEMIJI
D0 1505 L = I.N

M = MEMBER(L.JI
O I GO TO 1503

IF( M 0LT.
C
C POSITIVE INCIDENCE
C.
DO 1‘50? I = 1.1
FOQCFII.J) = FORFF(I.JI - RMI1.I.MI*FPP(I.M)
I - PMI?.I.MI*FPPI2.MI
P - PMI39I.M)*FPPI3.M)
QM(I.I.MI*Fpp(4-9MI

1:30? FORCEII+1.JI = FOPCFII+3.JI ‘-
- PMI2.I.M)*FPPI"3.M)

RMI3.I.M)*FPP(6.MI

I

2
GO TO 1:053

C

C NEGATIVE INCIDENCF?
C

I

503 M =-M

DO 1904 I = 1.1

FORCEII.JI = EORCFII.JI u RM(1.I.MI*FNN(I.M)

l - RMI2.I.MI*FNN(2.MI

P - RMI3.I.MI*FNN(3.MI
1504 FORCE(I+3.JI = FORCF(I+3.JI - RM(1.I.M)*FNN(4.MI

I - RM(?.I.M)*FNN(G.MI

2 - RM(3.I.MI*FNN(6.MI

I ‘305 CONTINUE
I ‘306 CONTINUE

c ----------------------- r
C GROUND ACCFLFRATION C
r -------------- - -------- r

IF! QUAKE I I710. I700
171!) XTIMFT = TIMF. + OTIME
DO 1730 K = 1.3

C

C INDEXIKI TO BE DETERMINED SUCH THAT

C DATA T(K. INDEXIKI-II .LT. XTIME .LE. DATA T(K. INDFX(K)I
C

I720 IF( XTIMF. .LE. DATA T(K. INDEXIKII I GO TO 1730

INDEX(KI = INDEX(KI + I
GO TO I720
C
C GROUND ACC. IN EAST-UR—QOUTH COORD. SYSTEM INTFRPOLATFD
LINEARLY BY TWO SUITABLE CONSECUTIVE READINGS

INDF’X(KI-II I*
DATA A(K. INDEXU<II - DATA A(K. INDEXIKI-II I/
DATA T(K. INDEXIKII — DATA T(K. INDEX(KI-II I
T(K) + DATA A(K. INDEX(KI-II )*AMPLIFY

XTIMF' - DATA T(K.

73D T(K)
1
?

17‘30 T(KI

C
C
I

AAAA

 

 

 

105

C

C GROUND ACC. IN GLORAL COORD. SYSTEM

C
GARB = ROT(1I*T(1I + ROTI2I*T(3I
T(3I =-ROT(2I*T(1I + ROT(II*T(3I
T(II = GARE

C

C GDOUND VELo. ACCo. AND DIS.

C I
00 1760 K 1.1 ;

X0 VFLIK) = x0 VEL(K) + 0.s*I x0 ACC(K) + T(K) )*nT
x0 ACC1K) = T(K) r
1760 x0 DIS(K) = x0 DIS(K) I

+

I XO VEL(KI + 005*XO ACC(K)*DT I*DT ,
I700 CONTINUE P

 

I ”'i.
‘1' "‘1'" .

O
I
I
I
I
I
I
I
I
I
l
I
I
l
I
l
I
I
I
I
I
I
I
I
l
I
I
I
I
I
I
I
l
I
l
I
I
I
I
I
I
I
I
D

Tl_

c FREE-JOINT ACCELERATIONS AND VELOCITIES c '
c ------------------------------------------ c

00 1612 J = 1.JFRFF 5 no 1612 I = 1.6 s GARR = 0

on 1611 K = 1.6
1611 GARE = GARB + RERRDNII.K.JI*

I I FORCE(K.JI + GTLOADCK.JI + REXT(K.JI I
IF( I .LT. 4 I GARE = GARR - XO ACCIII

VEL(I.JI = VELII.JI + O.S*( ACC(I.JI + GARE )*DT
ACC(1.JI GARR

1612 DIS(I.JI
GO TO ?000
END

DIS(I.JI + DDISII.JI

susROUTINF RMATRIY( J. ALPHA I
IDIMENSICWJ’WWAQQIIO). FQUW(3.1I. ~JUW(1.3I. ICCX(3.1I
DIMENSION HI3.3I. HJ(1.3). HJHI3.3I
RFAL MMAss. JMAss. JIM
COMMON RFsp0N16.6.s0;
COMMON/ TYPF/ DENSITY. E. G
COMMON/ RM/ RMI3.3.100)
COMMON/ DEADM/ DEADM(100I
COMMON/ COMP/ COMRI3.100>
COMMON/ LENGTH/ LENGTHIIDD) s REAL LENGTH
COMMON/ AREA/ ARFAX(IOO). AREAY(100). AREAZIInnI.
1 IXXIIOO). IYVI100). 172(100) s REAL IXX. IYY. I72
COMMON/ MFM/ NMEMIsn). MEMRFRIID.sn)
c
C wx. wY. AND wz
c

JOINT COORDINATFS OF JOINT MASS CFNTFR

JMASK§I= VXM = VWMN:= VZM = (I
IN = NMEM¢UH

106

D0 15 I = ION

M = MFMEEQIIQJI

TEMP : SIGNI Oo?E*ALDHA¢ M I

M = IABSI M I

MMASS(I) = 005*( DENSITY*AREAXIMI+DEADMIM) I*LENGTHIMI

VXM = VXM + TEMP§COMDII¢M)*MMASSIII
VYM = VYM + TEMP*COMPI2.MI*MMASSIII
sz = v7M + TEMP¥COMD(3.MI*MMASS(II
15 JMASS = JMASS + MMASRII)

JMASS : 1.0/JMASS

wx 2 VXM*JMASS

wv = VYMfiJMASS

WZ = VZM*JMASS

C
C EIM a BRANCH INERTIA MATRIX. JIM = JOINT INERTIA MATRIX
C

 

DO 20 L = 193 $ DO 20 K = 1.3
?O JIM(KQLI = 0

DO 30 I 3 ION

M = MEMEEQIIQJI

TEMP = SIGNI 0025*ALDHA0 M I

M = IAESI M I

VX = WX — TEMp*COMp(IOMI
VY 3 WY - TEMP*COMPI?9MI
VZ = W2 ~ TEMP*COMP(30MI
UX = RMIIQIQMI*VX + QMIIQZQMI*VY + RMII¢30M)*VZ
UY = RMIchoMI*VX + QMIZOEOMI*VY + PMI2039MI*VZ
UZ : RMI3OIOMI*VX + QMI3020MI*VY + RMI303OMI*VZ

BIM(1.2) = BIM(2.1) =-MMASS(I)*UX*UY

BIM(2.3) = BIM13o2) =-MMASS(I)*UY*UZ

BIMIBol) = BIM11.3) =—MMASSIII*UZ*UX

TEMP = 0.5*DENSITY*LENGTH(M)

BIM13.3) = TEMP*122(M)

BIMI2.2) TEMP*IYY(MI

BIMIIOII BIMI391I + RIMIPO?)

TFMP a ALPHA**?*LFNGTH(M)**1*

I I DENSITY¥ARFAXIMI+DFADMIMI 1/06.o

Rtmtsoa) BIM(3.3) + MMASS(I)*( UX**?+UY**2 ) + TFMP
BIM(2.?) BIM(2.2) + MMASS(I)*( UZ**?+UX**2 ) + TFMP
BIM(1.1) BIM<1¢1I + MMASS(I)*( UY**2+UZ**? )

DO 30 K = 1.3

Do 25 L = 1.3

25 XXXIKQLI = RMIIOKOMI*EIMIIOLI
I + PMI2OK¢MI*EIM(20LI
2 + RMIBQKQM)*EIM(39LI
DO 30 L = 103

30 JIM(KQLI = JIM(KOLI + XXXIK91I*PM(IOLQMI
I + XXXIK92I*RM(29LQM)
2 + XXXIK03I*RM(39L0MI

man

40

50

60
65

80

JIM

107

= INVERSE OF JOINT INERTIA MATRIX

CALL INVERSEI JIM93030HIIQIIQIoOE-QQJRQHII02)0HI193)

RESPON = JOINT RESPONSIVENESS MATRIX

HIIOII
HI203I
HISOII
HIIO?I
DD 50 K
DO 40 L

HI202I = HI303I = O
WX s HI392I =-WX
WY s HIIOBI B-WY
W2 5 HIPQII =-WZ
I03
103 S HJIKOLI = O S DO 40 I z 193
HJIKQLI HJIKQLI + HIKoII*JIMIIoLI
DO 50 L 103 S HJHIKQLI = 0 $ DO SO I = 1'3
HJH(K9L) = HJHIKOLI + HJIKIII*HII9LI
DO 65 K = 103
DO 60 L = 193
RESPONIKcLoJI =~HJHIK9LI
RESPONIK9L+30JI a RESPONIL+30K9JI =-HJIK9LI
RESPONIK+30L+39JI = JIM(K0LI
RESPONIKoKoJI : RESPONIKQKQJI + JMASS
PRINT 80. J. ALPHA. IIRESPONIK.L.J). K=1.6). L=1.6)
FORMATI*ORESPONSIVENESS MATRIX FOR JOINT*I3
* AND ALPHA =*F50?0 6I/EX6EI708II

END

SUBROUTINE FMATRIXI M I

C ..................... _- ........ C
GEOMETRICALLY LINEAR VERSION C
c--— ---------------- --—-——- ----- r

C

OI‘CIOIW

OI}

I

XII

XIZI
XI3I

HII

DIMENSION XI3I9 HI3I0 FLEXI893I9 FLEXIEIBI

COMMON/ AREA/ AREAXIIOOI. AREAYIIOOIo APEAZIIOOIQ
IXXIIOOI. IYYIIOOIo I7ZIIOOI S REAL IXXo IYY.

COMMON/ PLASTIC/ PLASTICI2OIOOI S LOGICAL PLASTIC

COMMON/ STIFF/ SIBOIOOIO BIZOIOOI

COMMON/ TYPE/ DENSITY. E. G

LENGTH OF POSITIVE~END PLASTIC PORTION
LENGTH OF NEGATIVEfiEND PLASTIC PORTION
LENGTH OF ELASTIC PORTION

I

XIII = XI?) = O

IFI PLASTICIIQMI I XIII PoO*BIIoMI
IFI PLASTICIZQMI I XI?) 200*EI19MI
XIBI = EIIoMI + EI?9MI - XIII - XI?)

I 2 DIS. BETWEEN POSo~END YIELD HINGE AND POS- END

I

Ifi

 

IZZ

108

C HIE) = DIS. RFTWFEN NIFGo-FND YIFLD HINGF AFID DOS. FND
C
HIII = OoE*XIII
HI?) = RIIOMI + RI70MI - Ooq*XI?I
HI3I = XIII
C
C FLEX. MATRICES FOR PLASTIC PORTIONS
C
DO 402 I = ID?
FLEXI79II = FLFXIEoII = O
TEMP = XIII/E
FLEXIIQII = TEMP/AREAXIMI
FLEXIEOII = TEMP/IVYIMI
FLFXI69II = TFMP/I7ZIMI
TEMP = XIII/G
FLEXIEOII = TEMP/ARFAYIMI
FLFXI30II = TEMP/ARFAZIMI
40? FLEXIAvII = TEMP/IXXIMI
C
C FLEX. MATRIX FOR ELASTIC PORTION
C
TEMP = XI3I/E
FLEXIIQQI : TFMP/ARFAXIMI
FLEXISQ?) = TFMP/IYYIMI
FLFXI603) = TEMP/I77IMI
FLFXI403I = XI3I/I G*IXXIMI I
FLEXI2¢3I = XI3I**?*FLEXI603I/300
FLEXI303I = XI3I**?*FLEXIEQ3I/Bon
FLEXI7v3I =-O.E*XI3I*FLEXI693I
FLEXIEQBI = OoE*XI3I*FLEXI593I
C
C FLEX. MATRICES REFERRED TO POSITIVE END
C
DO 404 J = 193
FLFXIQGJI = FLFXI?9JI
1 + HIJ)*( H(J)*FLFXI69JI — ?.D*FLFXI7.J)
FLFXIB-J) = FLFXI3¢JI
I + HIJI*I HIJI*FLEXI59JI + 700*FLFXIEOJI
FLEXI79J) = FLFXI70JI - HIJI*FLFXI69JI
404 FLFXIEQJI = FLFXIEQJI + HIJI*FLFXIE¢JI
C
C FLEXIEILITY MATRIX
C
DO 406 I = 10E
406 FLFXIBIII = FLEXIIQII + FLEXIIQ?) + FLFXIIQWI
C
C STIFFNESS MATRIX
C
SIIoMI = IoO/FLEXIEIII
SIAOMI = IoO/FLEXIEIAI

I

I

 

109

TFMD = FLFXIEI?I*FLEXIEI6I - FLEXIEI7I**?
SIPcMI = FLLXIEI6I/TFMP

SI69MI - FLEXIEI?I/TFMP

SITQMI =-FLEXIEI7I/TFMP

TEMP = FLEXIBI3I*FLEXIEIEI - FLEXIEIEI**2
SIOQM) = FLEXIEISI/TFMP

SISQMI - FLFXIEIBI/TFMP

SIEqMI =-FLFXIEIEI/TFMP

END

SUEROUTINE FMATRIXI M I

c ----------------------------------- c
c GEOMETRICALLY NONLINEAR VREGION c
O ----------------------------------- c

COMMON/ TYPE/ DENSITVQ F. G
COMMON/ STIFF/ SIEQIOOIQ RIPQIOOI
COMMON/ AREA/ AREAXIIOOIQ AREAYIIOO). AREAZIIOOI.

I IXXIIOOI. IYYIIOOIO IZZIIOnI $ REAL IXX. IYY.
C
C CONSTANT FACTORS OF FLFMFNTS OF ELASTIC STIFFNESS MATRICES
C
X = EIIoM) + EIPOMI $ TEMP = F/X
S(SoM) = TEMP*IYYIM) $ S(EoMI = SISoMI/X
SI6¢MI = TEMP*IZZIM) s 9(7QMI = SI6oMI/X
S(loMI = TEMP*ARFAXIMI
SIZoMI = 200*SI79MI/X
SI3QM) = ?oO*SI80MI/X
SIaoM) = G*IXXIMI/X
END
SURROUTINE SMATRIXI EFTA. MLINFAR. GLINFAR )
c -------------------------------------------- — ---------- c
C GEOMETRICALLY LINEAR VERSION. PARAB. YIELD FUNCTION C
c ------------------------------------------------------- c

LOGICAL MLILEAR. GLINEAR. ITERATE

DIMENSION DWI6). FRIG). VI6.2I. HIE)

DIMENSION RHII2.1OO). CHIIPI. 88(2). EEI2.2I. XNI6)
DIMENSION SSI606I9 GGIOQEIQ SGI602I¢ XXI296IQ YYI696I
DATA (SS = 36(OII. (66 = 12(0))

COMMON/ TIME/ TIME. DT

COMMON/ RM/ RMI3.3.1nn)

COMMON/ SLIMIT/ :LIMITI6.IOO)

COMMON/ JDJN/ JPIInfi). JNIIOOI

COMMON/ STIFF/ S(R.IOOI. HIR.1OO>

 

‘_I “I“

110

SI7F/ MEMBERS. JOINTS. JEREE
FORCE/ FPPIéoIOOI! FNNIfiOIOOI
LENGTH/ LENGTHIIOOI 3 REAL LENGTH
PLASTIC/ PLASTICI2OIOOI $ LOGICAL ELASTIC
COMMON/ YIELDED/ YIELDEDIEQIOOI S LOGICAL YIELDED
COMMON/ DIS/ DISI6¢EOI9 DDISIéoEOIo
1 PDISPIfiolOO). PDISNI691OO)
DO EEO M: IQMEMEFDQ
c ------------------- ~ ------------- c
INCREMENTAL MEMBER DISTOPTION
C ------------- ‘ ------------------- C
JPIMI
JNM JN(M)
DO 305 I
DWII)
DO 305 J
XNII+3I =
DWII+3I =
I + RMII.J0M)*I DDISIJ+30JPMI~DDISIJ+39JNM) )
DWII) DWIII + RMII.J9MI*I ”DISIJOJPM)-DDISIJQJNM) I
DWI?) DWI?) + LENGTHIMI*XNI6I
DWI3) DWIqI m LENGTHIMI*XNIEI
c ------------------------------------------- O
C POGO-FND IF MFMRFD M NOT VIFLDING

COMMON/
COMMON/
COMMON/
COMMON/

103
DWII+3I
I03
XNII+3I
DWII+3I

XNII+3I O

+ RMIIOJQMI*DDISIJ+39JNMI

10G

COR. PLASTICI2.M) I
OANDO pLASTICIszI I
GO TO 440

PLASTICIIQM)
PLASTICIIoMI
N oGTo O I

IEI
IEI
IFI

N

FPIII
FRI2I
FPI3)
FPIAI
ERIE)

FPPIIQM)
FPPIzoM)
FDPI3QMI
FDP(4.M)
FPPISQM,

GIIoMI*Dwt1)
S(PoM)*DW(?I
§I3.M)*DWI3I
GIA.M)*DMIAI
SIQQMI*DWIQI

+
+

+

SI70M)*DWI6)
SIBQMI*DWI5I

SIEQM)*DWI1)

+.+.+4.+ +

EPPIficM) + SI79MI*DWI?I
a 19?
AESIEPI1))*SLIMITIIOM) + IFPI4I*SLIMITI49MII**?

I + ARSIEPIEI + EIJQMI*FPI3))*SLIMITI59MI

+ ABSIFPIéI RIJqMI*FpI?II*SLIMITI60MI
IEI MLINEAP I EDD. 7OO
CONTINUE
C --------------------------------------- C
C ano~FND FORCE IF MEMEED M YIELDING

c --------------------------------------- c

FPIfiI SIAQMI*OWI6I
DO A?O

CHIIJ)

J

420

1.0

EI?0MI
oAND. PLASTIFI29M) I

EIloMI G HI?)
ONOTQPLASTICIIOMI

HIII HIP)

YIELD SURFACE NORMAL:

111

DO 50? J = ION

VIIOJ) = SIGNI SLIMITII9MIO FPPIIQM) )

VI4OJI = 200*FPPI40MI*SLIMITI49MI**?

VIEOJ) = SIGNI SLIMITISOMIO EPPISOMI+HIJI*EPPI39M) )
50? VI6OJI = SIGNI SLIMITI6OMIO FPPI69MI-HIJI*FPPI2OMI I
C
C STIFFNESS MATRIX
C

DO 503 I = 196
503 SSIIOII = SIIOM)

SSI692I 55(296I = SI79MI
SSIR.3I - SSISOEI - SIEOMI

GG MATRIX

'J'IDDO

OE ITERATE a oEALSE.
DO 510 J = ION
GGIIOJI = VIIOJI
GGIZOJ) =-VI6OJI*HIJ)
GGISOJ) = VISOJ)*HIJI
DO 510 I = 406

E10 GGIIOJ) = VIIOJI

C

C SG = SS*GG

C
DO 511 I a 196 $ DO E11 J = ION S SGIIQJI = D
DO 51! K = IO6

511 SGIIOJ) 8 SGIIOJI + SSIIOKI*GGIKOJI

C EE : INVERSE OF TIGGI¥SS¥GG T MEANS TRANSPOSE
DO 512 I 8 ION 5 DO E1? J = ION 5 FEIIOJ) = 0
DO 512 K = 196

512 EEIIOJ) = EEIIoJI + GGIK9II*SGIK9JI
CALL INVERSEI EEO N9 N9 XNIII9 IoOE-QO L1 XNI3I9 XNIEI I

C
C XX 8 EE*TISG) = EE*TIGGI*TISSI
C

DO E14 I 8 ION 5 DO :14 J a 196 5 XXIIOJI = 0
DO 514 K s ION

514 XXIIoJ) = XXIIOJ) + FEIIoKI*SGIJ9KI

C

C YY = SG*XX = SS*GG*EE*TIGGI*TISSI

C
DO 518 I = IO6 5 DO 518 J 3 1'6 5 YYIIoJI = 0

DO 518 K a ION
51R YYIIOJI = YYII9J) + SGIIQKI*XXIK9JI
C
C POSITIVE-END FORCE
C

112

DO :3? I = 196
FP(II = O
DO 551 J = 196
551 FPIII = FDIII + ( SS(IOJI-YY(IOJI I*OW(J)
55? FP(II = FPIII + FPDCIOMI
DO 553 J = IO?
CHI(JI = ABSIFP(1))*SLIMIT(19MI + (FP(4I*SLIMIT(4OM))**2
I + ABSIFPISI + B(J9MI*FP(3II*SLIMIT(59M)
P + ABSIFPIGI - B(J9MI*Fp(2II*SLIMIT(69MI ‘ 100
IF( PLASTIC(JOMI oANDo I CHIIJI oGTo O I I ITFRATF = oTRUFo
553 CONTINUF
C
C YIELD SURFACE NORMALS
C
IF! ITEPATF I 5709 GOO
570 L009 = LOOP + I
DO 572 J = ION

 

VllOJI = SIGN( SL1MITI19MIO FP(1I I
V(4OJI = 200*FP(4I*SLIMITI49MI**2
V(59J) SIGN( SLIMIT(59MI9 FP(GI+H(J)*FP(3) I

57? V(6OJI SIGN( SLIMITI69MIO FP(6I-H(J)*FP(?I I
IF( LOOP 0LT. 4 I GO TO 505
PRINT E819 TIME s CALL FXIT

581 FORMAT(*-LOOP IS TOO LARGE AT TIME =*F8.5I

590 CONTINUF

c -------------------- — ------ c
C TEST FOP ELASTIC QETURN C
C—-— ------------------------ C
DO 640 J 8 ION
BB(J) = 0

DO 620 K = IO6
6?O BB(JI = BB(JI + XXIJOKI*DW(KI
IFI BBIJI OGT. O I GO To 640
IF( oNOTopLASTICIIOMI OAND. pLASTIC(29MI I J
FLASTICIJOMI 8 OPALSEO
IFI BETA oGTo O I CALL FMATRIXI M I
L 8 M $ IFI J CFO. ? I L 3 -M
pQINT 6?49 LO TIMEO CHIIJI
624 FORMAT(*-MEMBER*I3* PFTURNS TIMF=*FBOQ* CHI=*F1?OBI
640 CONTINUE
C
C pLASTIC MEMBER DISTOQTION
C

II
\I

DO 660 I = IO6
DO 660 J : ION
660 POISPIIOMI = pDISP(IOMI + GG(IoJI*RR(JI

c--- ------------------ c
c TFQT FOR YIELDING c
C—-- ---------- ~ ----- --n

700 DO 790 J 3 192

113

IF! PLASTIC(J9M) ) GO TO 780
IF! CHIIJ) .LTo o I GO TO 760
PLASTIC(J9M) = YIFLDFD(19M) = .TRUF.
IF( BETA .GT. 0 ) CALL FMATRIX( M )
L = M 5 IF( J .FQ. P I L = -M
PQINT 7?4. L9 TIMF. PHIIJOM)
7pa FORMAT(*-MEMBEP*13* YIELDS TIME=*FR.R* PHI=*F1?98I
GO TO ann
760 RHIIJOMI = CHIIJ)
van CONTINUF

c--- ------------------ c
c MFMBER END FODCES c
c——— —————————————————— c

BOO FPR(R9M) = IPIR)
FPP(69M) = FP(6)
FNN(89M) =-FP(=) - LFNGTH(M)*FP(3)
FNN(69M) =-FP(6) + LFNGTH(M)*FP(2I
DO 830 I = 194
FPPIIOM) = FP(II

810 FNN(IOM) =~FPIII

BRO CONTINUF
FND

SURROUTINF SMATRIX( RFTAO MLINEAR9 CLINFAR I

c -------------------------------- — ------------------------- c
C GFOMETRICALLY NONLINEAR VRESIONO ELLIP. YIELD FUNCTION c
c ---------------------------------------------------------- c

LOGICAL ENDPO ENDNO ENDSO MLINEAR9 GLINEARO ITERATE
DIMENSION DUP(6I9 DUNI6IO FPIfiIO FNI6IO VpIéIO VN(6I
DIMENSION PHI(?O1OO)9 CHI(?I9 BB(ZIO FE(?O?I9 XNIGI
DIMENSION GRI692IO GNIéOPIO FGRT(?96IO FGNTIPOéI
DIMENSION GPEGRTI696IO GDFGNTI696IO GNFGNTI696)
DIMENSION SRPC606I9 QRNHSOFIIO RNNIGOEI

DATA (9RD = 36(OIIO (QPN = 16(OIIO (SNN = 36(0II
COMMON/ TIME/ TIME. DT

COMMON/ RM/ RM(303OIOOI

COMMON/ SLIMIT/ SLIMITIéOIODI

COMMON/ JPJN/ JPIIOO)9 JN(IOOI

COMMON/ STIFF/ SIBOIOOIO B(?OIOOI

COMMON/ SIZE/ MEMBERSO JOINTSO JFREE

COMMON/ FORCE/ FRRI6OIOOIO FNNI6OIOOI

COMMON/ LENGTH/ LENGTHIIOOI 5 REAL LENGTH

COMMON/ PLASTIC/ RLARTICIEOIOOI $ LOGICAL pLASTIC
COMMON/ YIELDED/ YIELDFD(?OIOOI $ LOGICAL YIELDEO
COMMON/ DIS/ DIQISOGOIO DDIQI6OQOIO

I RDI$P(6OIOOIO RDISNIéOIOOI

DO 880 M I IOMEMRERS

 

114

r ---------------------------------------- r
C INCREMENTAL MEMBER END DISPLACEMENTS C
C ———————————————————————————————————————— r

JPM = JPIMI
JNM 2 JN(M) .
DO I03 I = 193

DUPIII = DURII+31 = DUN(II = DUNCI+3I = O
DO 10; J 3 193
DUP(II = DUPIII

+
DURII+3I = DUPII+1I + RMIIOJOMI*DDI§(J+3OJDMI
DUNIII = DUNIII + RMIIOJOMI*ODIS(J9JNMI

IO; DUNII+1I = DUNII+3I + RMIIOJOMI*DDI§(J+39JNMI

r —————————————————————————————————————— — -------- c

C FACTORS NEEDED IN MEMBER STIFFNESS MATRICES C

C--— ----- - --------------------- - ------ ----------c a
IFI GLINEAR I GO TO R60
9 = FNNIIOMI
ABSR = ARSI P I

RMIIodoMI*DDIR(J9JPMI
,‘é

 

‘Il

 

C
C TFMP = 1 PEP CENT OF fiMALLE9T FULER LOAD
C
TFMP = n.n?€*$(8oMI
IF( ARSP .LT. TFMD I GO TO 260
PHI? = 900TI ARQP/QI79MI I
PHIB = SQPT( ARSP/S(BOMI I
IF( P I 2209 2609 240
2PO C? = COSt Ple I s S? = SINI DHI? I
C3 = COS! PHIB I s Si n SINI PHI3 I
GO To 280
240 C? 0.=*EXP( PHI? I s 57 3 EXPI—PHI? I
C1 = OoG*FXP( PHIfi I $ 91 a EXP(-PH13 I
C? = C? + O.G*S? g Q? = C? - 5?
C1 = C3 + 005*93 s §1 = C1 - $3
?=0 TFMD = -SIGN( 1.09 D I
CARE = TEMP/( ?.0*( I.O~C? I/PHI? - TFMD*§P I

H?! = GADB*( S? — C?*DHI? I
HP? = GAQR*I PHI? - S? I
H23 = HPI + H2?

GARE = TFMp/( PoO*( 1.0-C3 I/pHIB - TFMR*83 I

H31 = GAQB*( $3 - C3*PH13 I
H32 = GARB*( PHI3 - 93 I
H33 = H31 + H3?

RH02 = I RMI3OIOMI*( DIS(IOJRMI-DIQIIOJNMI

I

I + RM(39POMI*( DIS(HOJPMI-DI§(29JNMI I
2 + RMI3O19MI*I DIRI3OJPMI-DIRI39JNMI I I/LFNCTH(MI

pH03 = ( RM(?91OMI*( OIR(19JNMI-DIQ(19JRM) )

I + DM(?9?OMI*( DI:I?QJNMI“DIQI?QJDMI I
P + RM(293OMI*( DIS(fiOJNMI-DIRI3OJRMI I I/LFNGTH(MI

TEMP = P/LENGTHIMI
GO TO PRO

115

P6O CONTINUE
H31 = HP] = 4.0 5 HR? 3 HP? = 290 5 H31 = H23 : 6.0
RHOS = DHO? = TEMP = O
?RO CONTINUE
C-- --------------------------- C
C MEMBER STIFFNESS MATRICES C
C ----------------------------- C
QPPIIOII = QNNIIOII = QIIOMI
QDPIAOAI = QNNIAOAI : RIAOMI
gPPIROQI = RNN(QQGI = H71*9(ngI

 

RPP(696I = SNN(696I H91*§(69MI
SPP(1.2I : SPP(?OII = SNN(IO?I = SNN(?OII = RHO1*SPP(191I
SPPIIO3I = SPPI3O1I = SNNI193I = SNN1391I =~RHO?*SPP(IOII
SPPIZOBI = SPPIBOZI = SNN(293I = SNN(392I = RH03*SPP(193I
SPPC292) = SNNIEO?) = TEMP + RH03*SPP(1921 + H23*5(2OMI
SPRI3O3I = SNNI303I = TEMP - RH02*SPP(193I + H33*S(39MI
DO 310 J = 193 5 DO 110 I = IO3
110 SDNIIOJI =~§PPII9JI
RDNIQO4I g—QPP(4.4)
SPNIQORI = H3?*9(:9M) $ TEMP = H?3*Q(7OMI
RPNI606I = H??*§(69MI $ CARR = H11*<(POMI
9PP(?96) = SPPIGO?I = quIPC6I = TEMP
SNN‘296I = SNNI69?I = SPNISOPI =-TFMP
SDRI3OSI = SPPI593I = SPNISOSI =-GARB
SNN(39€I = SNNI593) = SPNISOEI = GARE
c ---------------------- — ---------------- c
C END FORCES IF MEMBER M NOT YIELDING C
C ----------------------------------- ----C
400 CONTINUE
ENDP = EmN = ENDR 2 oFALSEo
IFI PLASTIC(IOMI oANDo PLASTIC(?OMI I FNDS = oTRUFo
IEI oNOToFNDS oANDo PLASTIC(19MI I ENDP = oTRUFo
IEI oNOToFNOS oANDo DLASTICIZOMI I ENDN = oTRUF.
IFI ENDQ oORO ENDP oORo ENDN I GO TO 480
D0 440 I = 196
FRIII = FNIII = O
DO 4?O J = IO6
FP(II = FP(II + SPPII.J>*DUP(J) + QPN(IOJI*DUN(JI
4?O FN(II = FN(II + SPNIJOII*DUP(JI + SNN(IOJI*DUN(JI
FRIII : FD(II + FPPIIOMI
440 FNIII : FN(II + FNNIIOMI
CHIIII = (FPIII*§L1M1TI19MII**? + (FPI4I*SLIMIT(4OMII**?
I ‘ 1.0 + (FPIQI*§LIMITIQOMII**? + (FPI6I*§LIMIT(6OMII**7
CHIIp) = (FN(1)*SLIMIT(19MII**? + (FNI4)*SLIMIT(4.M))**?
I - 1.0 + IFN‘RI*QLIMIT(ROMII**? + (EN(6I*9LIMIT(69MII**?
IEI MLINFAR I BOO. 7nn
49O CONTINUE
C ----------------------------------- C
C END FORCES IF MEMBER M YIELDING C

C-—- ---------------------------- -‘--C

C
C
C

SO]

O

C
C

60:

510

613

5?0

6?3
813

0(WCIO

116

LOOP : 1
YIELD SURFACE NORMALR

IF( FNDN ) GO 70 :n1
VP(1I = Fpp(1.MI*9LIMIT(19MI**P

VP<4I = FPP(4OMI*SLIMIT(AoMI**7
VP(SI = FDP<5.MI*<LIMIT¢a.MI**?
the) = FPPI69MI*SLIMIT(69MI**2

IF( ENDP I GO TO R06

VN(II = FNN(I.M)*SLIMIT(19MI**2
VN(AI = FNNI4.MI*SLIMIT(4OMI**?
VN(EI = FNNISOMI*SLIMIT(R.MI**2
VN(fiI = FNNI69MI*RLIMIT(69M)**?

GPO GNO AND FE MATRICFQ

ITFRATE : QFALQE.
IF( ENDN I GO TO :11
DO RIO J = 196

GDIJqu = SPPIJ.1I*VD(1I + <PR(J94I*VP(41
l + SPPIJ98I*VDIRI + SODIJO6I*Vp(6I
GN(J91I = SPN(19JI*VP(1I + SPNtanI*VP(aI
1 + SPN(RoJI*VP(5I + SPN(6.JI*VD(6I
EF¢1~II = GP(1¢1)*VP(1I + GP(491)*VP(4)
1 + GP(=.1I*VP(RI + GP(691I*VP(6I
IF( FNDP I GO TO =23
DO GPO J = 1.6
GP(J9PI = SPNIJ91I*VN(1I + QPNIJ94I*VN(4I
I + SPNIJ9=I*VM(GI + GPN(J96I*VN(6I
GNIJOPI = SNNtJ.1)*VN(II + 9NN¢J94I*vN(4I
1 + SNN(J98I*VM(RI + QNN!J.6I*VN(6I
EF(?92I = GN(1o2I*VN(1I + GN(4.2I*VN(AI
1 + GNIR.?I*VN(5I + GNI6.2I*VN(6I

IFI FNDN I GO TO :33

55(19?) GP(19?I*VD(1I + GD(49?I*VP(aI
1 GP(=.2I*VPIGI + GPI69?I*VP(6I
EF(291I - FE<19?I

CALL INVERSF( PE. 2. 2. XNIII. 1.6F-9. L.

|+|l

IF( FNDD I EF(191I = Ion/FF(191I
[Ft ENDN I EF(?9?I = 1.0/FF(29?I
GDEGDT GP4EE§TIGPI T MEANQ TRANqDOSF

GDEGNT = GP*FF*T(GNI GNFGNT = GN*F=*T(aN)

IFI ENDS I GO TO Ran

J = 1 $ IF( ENDN I J = P
DO 535 K = 196

EGPTIJOKI EEIJOJI*GP(KOJI
EGNT(JOKI = EEIJoJI*GN(K9J)

XN(3I9

XN(:)

I

£34 a .Auau "k"

_jld"""' "

 

117

DO 83: I = IO5
GPFGPTIIOKI
GPFGNTIIOKI
Q13 GNEGNTIIOKI
GO TO Gan
“I40 00 PAR K = 1‘6
EGPTIIOKI EEIIOII*GP(KOII
EGPT(29KI EEI2O1I*GP(KOII
EGNTIIOKI EEIIOII*GN(KOII
EGNTIQOKI EFI?OII*GN(KOII
DO E45 I = IO6
GPEGPTIIOKI
GPFGNTIIOKI
54R GNEGNT(IOKI

GPIIOJI*EGPTIJ9KI
GPIIQJI*FGNTIJOKI
GNIIOJI*FGNT(JOKI

FEIIO?I*GPIK9?I
FE(?O?I*GP(KO?I
FEIIOPI*GN(KO?I
FEI?O?I*GN(K9?I

IIIIH
++++

GP(191)*EGPT(IOKI + GPIIO2I*FGPT(?OKI
GP(191I*FGNT(IOKI + GP(192I*EGNT(29KI
GN(191I*EGNT(IOKI + GN(IO2I*EGNT(?OKI

C

C END FORCE:

C

SRO DO 5‘2 I = IO6

FPIII = FN(II = O
DO 551 J = 196
FP(II FPIII
I

551 FNIII
I

SPPIIOJI-GPEGPTIIOJI I*DUP(JI
SPN(IOJI-GPEGNT(IOJI I*DUN(JI
SPN(JOII-GPEGNT(JOII I*DUP(JI
SNNIIOJI-GNEGNT(IOJI I*DUN(JI

AAA“

+
+
FNIII +
+
+

FDIII = FDIII FDPIIoMI

SR? FN(II = FNIII + FNNIIOMI
CHI(1I = (FPI1I*SLIMIT(19MII**P + (FP(4I*SLIMIT(49MII**P
l — 1.0 + (FP(8I*SLIMITIG.MII**? + (FPI6I*§LIMIT(6~MII**2
CHIIpI = (FNI1)*SLIMIT(1.MII**2 + (FN(4I*SLIMIT(4.MI)**?

I - 1.0 + (FN(€I*SLIMIT(S.MII**? + (FN(6I*SLIMIT(6OMII**?
IF( PLASTICIIoMI .AND. (CHI(1I .GT. 0) I ITERATE = .TRUE.
IFI PLASTIC¢29MI oANDo (CHI(2I oGT. OI I ITERATE = oTRUE.
C
C YIELD SURFACF NORMAL:
C
IF( ITFPATE I 5709 890
570 LOOP = LOOP + 1
IF( FNDN I GO T0 ‘71
VPIII = FPIII*SLIM1T(1oMI**7
VP(4I = FP(4I*§LIMIT(49MI**?
VPIOI = FDI5I*€LIMIT(RoMI**P
VP(6I = FP(6I*9LIMIT(69MI**2
IF( ENDP I GO TO R7?
571 VN(II = FN(1I*SLIM1T(19MI**?
VN(AI FN(4I*SLIMIT(40MI**?
VN(SI FN(5I*§LIMIT(%oMI**2
VN(6I = FNI6I¥SLIMIT¢69MI**P
572 IF! LOOP .LT- 4 I Go TO 805
PRINT 5819 TIMF % CALL EXIT
58! FORMAT(*-LOOP IS TOO LARGE AT TIMF :*FR.OI

 

118

600 CONTINUE

c --------------------------- c
c TEST FOP ELASTIC pETupN c
c --------------------------- c

DO 640 J = IO2
IF( PLASTICIJOMI I 6109 640
610 PBIJI = O .
DO 620 K = IO6
6?O BB(J) = RRIJI + EGPTIJOKI*OUP(KI + EGNTIJOK)*DUNIKI l -
IF( BB(JI .GT. n I GO TO 640 ”‘1
pLASTICIJOMI = OFALgFO fi
L=M$IFIJOFOO?IL="M f
PRINT 6?49 L9 TIMEO CHIIJI
624 EORMAT(*—MEMBER*13* RETURNS TIME=*F896* CHI=*E1298I
640 CONTINUE
C
C PLASTIC MEMBER END DISPLACEMENTS IN GLOBAL COORD9 SYSTEM
C

 

I'l'I-"u‘ (I! 9 .“ -. '

IEI ENDN I 6609 66O
66O DO 666 I = 193
PDISPIIOMI
DO 666 J = 193
666 PDISPII+39MI = PDI6PII+19MI + RM(J9IOMI*PR(II*VDIJ+3I
IEI ENDP I 67OO 66D
66D DO 666 I = 193
PDISNIIOMI = PDISNIIOMI + RMIIOIOM)*RR(2I*VN(II
DO 665 J = 103
666 PDISNII+39MI = PDISNII+3OMI + RMIJOIOMI*RR(2I*VNIJ+3I
67D CONTINUE

PDISPIIoMI + QM(1919M)*QR(1)*VD(1I

c ------------ ~ -------- c
c TEST FOP YIELDING c .
c --------------------- c

700 DO 780 J = 192
IF( PLASTICIJ9MI I GO TO 780
IF( CHI(JI 9LT. n I GO TO 76“
YIFLDEDIJomI : PLA¢TIC(J9MI : .Tnuc.
L = M 6 IEI J QED. 7 I L = ‘M
PRINT 7249 L9 TIMF9 PHIIJ9MI
794 FOPMAT(*~MEMBER*I1* YIELOS TIME=*FR.P* PHI=*F1?9RI
GO TO 400
760 PHIIJOMI = CHI(JI
780 ’ CONTINUE

c--— ------------------ c
c MEMBER END FORCES c
c--- ------------------ c
800 on 930 1 = 196
FDP¢I9MI = FP(II

81D ENNIIOMI
BRO CONTINUE
END

FNIII

"I711111171111"1111“