S-Frame Theory Manual

S-Frame
Theory
Manual
Table of Contents Theory Manual
Table of Contents
Introduction . . . . . . . . . . . . . . . 6
Nomenclature . . . . . . . . . . . . . . . 8
Analysis procedures . . . . . . . . . . . . . 12
Linear static analysis . . . . . . . . . . . . . 12
Analysis description . . . . . . . . . . . . . 13
Implementation details . . . . . . . . . . . . 14
Free-vibration analysis . . . . . . . . . . . . . 15
Stressed free-vibration analysis . . . . . . . . . . . 21
Linear dynamic analysis . . . . . . . . . . . . 23
Implementation details of Newmark algorithm . . . . . . . 25
Implementation details of adaptive time-step size algorithm . . . . 27
Selection of constant time-step size . . . . . . . . . . 28
Selection of minimum and maximum time-step sizes . . . . . . 29
P-delta analysis . . . . . . . . . . . . . . 29
Linear buckling analysis . . . . . . . . . . . . 31
Buckling - Stressed Free Vibration - P-delta . . . . . . . . 34
Response spectrum analysis . . . . . . . . . . . 35
Nonlinear static analysis . . . . . . . . . . . . 43
Newton-Raphson Method . . . . . . . . . . . . 45
Convergence Criteria . . . . . . . . . . . . . 45
Other procedures . . . . . . . . . . . . . . 48
Skyline storage scheme . . . . . . . . . . . . 48
Mesh optimization . . . . . . . . . . . . . 51
Application of boundary conditions . . . . . . . . . 52
Constraints (coupled degrees of freedom) . . . . . . . . 54
Damping effects . . . . . . . . . . . . . . 57
Checking for non-structural degrees of freedom. . . . . . . 59
Graph interpolation . . . . . . . . . . . . . 59
Lumped mass VS. consistent mass formulations . . . . . . . 60
Number of elements required . . . . . . . . . . . 61
Friday 28 April 2000 – 14:40

Element library . . . . . . . . . . . . . . 62
3D Beam . . . . . . . . . . . . . . . . 64
Stiffness, geometric stiffness, and mass matrices . . . . . . . 65
3D Thin Shell . . . . . . . . . . . . . . . 68
Stiffness and mass matrices . . . . . . . . . . . 68
Plate elements . . . . . . . . . . . . . . 68
Membrane elements . . . . . . . . . . . . . 68
Aspect ratio . . . . . . . . . . . . . . . 69
3D Truss . . . . . . . . . . . . . . . . 69
Stiffness and mass matrices . . . . . . . . . . . 69
Two-noded spring . . . . . . . . . . . . . . 70
Stiffness matrix . . . . . . . . . . . . . . 71
Tension-only members . . . . . . . . . . . . 71
Compression-only members . . . . . . . . . . . 71
Spring element . . . . . . . . . . . . . . 72
Stiffness matrix . . . . . . . . . . . . . . 72
Hook element . . . . . . . . . . . . . . . 72
Gap element . . . . . . . . . . . . . . . 73
Unit conversion factors . . . . . . . . . . . . 74
Memory requirements . . . . . . . . . . . . 78
Dynamic memory allocation requirements (RAM) . . . . . . 78
Storage requirements . . . . . . . . . . . . . 78
References . . . . . . . . . . . . . . . . 80
Friday 28 April 2000 – 14:40

Friday 28 April 2000 – 14:40

S-Frame Theory Manual Chapter 1 : Introduction : 6
1 Introduction
This document describes the theoretical and computational
aspects of the S-Frame finite-element analysis engine. It is hoped
that this information will enable the user to make full use of
S-Frame’s capabilities while being aware of its limitations.
D:\softek docs\book\Theory\Theory.csc.fm
7 : Chapter 1 : Introduction S-Frame Theory Manual
Friday 28 April 2000 – 14:43

S-Frame Theory Manual Chapter 2 : Nomenclature : 8
2 Nomenclature
Following is a list of the symbols used throughout this manual.
A cross-sectional area.
s
Ai shear area normal to direction i.
c generalized damping.
c cr critical generalized damping.
E Young’s modulus.
Fi Dominant modal base shear in ithdirection.
G shear modulus.
J jerk or torsional moment of inertia.
ˆ
J normalized jerk.
Jcr critical jerk.
k generalized stiffness.
L element length.
m generalized mass.
me element mass.
ˆ
m modal mass.
%M
j
percent of total mass participating in jth direction.
Pe axial load (positive for a tensile load).
Ri RSA Base Shear in ithdirection.
Si Scaling factor in ithdirection.

t time.
Ti Code base shear in ithdirection.
x generalized displacement.
wh half bandwidth.
α, β Rayleigh damping parameters.
α, δ Newmark parameters.
∆t time step.
9 : Chapter 2 : Nomenclature S-Frame Theory Manual
θ phase angle.
λ load factor or shift factor or scale factor.
λ cr critical load factor.
j
µ modal participation factor in jth direction.
ξ damping ratio ( c ⁄ ccr )
ω circular natural frequency.
ω shifted circular natural frequency.
Ω frequency of disturbing force.
[A] a square matrix.

[ c ]s system modal damping matrix.
[ C ]e element damping matrix.
[ C ]s system damping matrix.
ˆ
[D] a diagonal matrix.
[ k ]s system modal stiffness matrix.
[ K ]e element stiffness matrix.
[ K ]s system stiffness matrix.
[ Kr ] reduced system stiffness matrix.
)
[ K ]s a shifted system stiffness matrix.

[ K eq ] equivalent system stiffness matrix.
[ Kg ]e element geometric stiffness matrix.
[ Kg ]s system geometric stiffness matrix.
[ K go ]s system geometric stiffness matrix at reference
loading.
[L] a lower matrix.
[ m ]s system modal mass matrix.
[ M ]e element mass matrix.
[ M ]s system mass matrix.
[ Mr ] reduced system mass matrix.
Friday 28 April 2000 – 14:43

S-Frame Theory Manual Chapter 2 : Nomenclature : 10
ˆ
[ ]s a system matrix after the application of boundary
conditions.
{a} a vector.
{b} a vector.
{ F }e element internal load vector.
{ F }s system internal load vector.
j
{I} unit vector in jth direction expanded to system
dimension.
{J} system jerk vector.
{ R }e element applied equivalent load vector.
{ R }s system applied load vector.
{ R }s
i
system applied load vector for ith load
case/combination.
{ R eq } equivalent load vector.
{ W }s relative system displacement vector.
j max
{ w }i maximum relative displacement contribution of
mode i due to ground motion in the jth direction.
max
{ w }i maximum relative displacement contribution of
mode i due to ground motion.
max
{W} maximum relative displacement due to ground
motion
{ X }e element displacement vector.
{ X }s system displacement vector.
·
{ X }e element velocity vector.
·
{ X }s system velocity vector.
··
{ X }e element acceleration vector.
··
{ X }s system acceleration vector.
{ Z }s ground displacement vector.
·
{ Z }s ground velocity vector.
11 : Chapter 2 : Nomenclature S-Frame Theory Manual
··
{ Z }s ground acceleration vector.
{ ˆ }s a system vector after the application of boundary
conditions.
{ η }s generalized system displacement vector.

·
{ η }s generalized system velocity vector.
··
{ η }s generalized system acceleration vector.
[Λ] natural frequencies and eigenvalues matrix
(diagonal).
[Φ] mode-shapes or eigenvectors matrix.
{φ} mode shape or eigenvector or buckling mode.
{ψ} normalized mode shape or eigenvector.
[Ψ] normalized modal matrix or eigenvectors matrix.
numdir number of ground motion directions.

num1Delm number of 1D elements (beams, trusses, etc.)
num2Delm number of 2D elements (shells, plates, etc.)
numelm number of elements.
numbc number of boundary conditions (constrained
degrees of freedom).
numdof number of degrees of freedom in structure.
numitrvec number of iteration vectors.
numldcas number of load cases.
numldcmb number of load combinations.
numnatfrq number of natural frequencies.
numnod number of nodes.
numRSAldcasnumber of RSA load cases.
numtimstp number of time steps.
tstval test value.
RMS root mean square.
Friday 28 April 2000 – 14:43

S-Frame Theory Manual Chapter 3 : Analysis procedures : 12
3 Analysis procedures
This section briefly describes the theoretical basis of the various
analysis procedures available in S-Frame. S-Frame uses the
displacement method of the finite-element method. In this
method, the structure to be analyzed is approximated by an
assembly of structural regions (elements) connected at a finite
number of points (nodes) to ensure that the displacements are
continuous. Once the equilibrium equations for each element are
known
{ F } e = { R }e (1)
then their contribution to the behavior of the overall structure

can be accounted for by assembly of the element equations (Eq.
1) using standard matrix procedures
numelm
∑ ( { F }e = { R }e ) (2)
e=1
to yield the system equations
{ F }s = { R }s (3)
where the internal load vector, for linear analysis, is defined as
linear static analysis

{ F }s = [ K ]s { X }s (4)
linear dynamic analysis

·· ·
{ F } s = [ M ] s { X } s + [ C ] s { X } s + [ K ]s { X } s (5)
nonlinear static analysis

{ F }s – { R }s = { 0 } (6)
Linear static The purpose of the static analysis is to determine the

displacements and stresses due to time-independent loading
analysis
conditions under the following assumptions:
13 : Chapter 3 : Analysis procedures S-Frame Theory Manual
1. Stiffness effects and applied loads do not depend on time.
2. Inertial and damping effects are ignored.
3. Static acceleration fields, such as gravity, may be included

(i.e. dead loads).
4. Time independent loads, displacements, pressures and

temperature effects may be applied.
Analysis description The system equilibrium equations for the linear static analysis
are
[ K ] s { X }s = { R }s (7)
where the system load vector { R } s includes the contributions of:
1. Applied nodal loads.
2. Loading due to static acceleration fields (such as gravity).
3. Element thermal/pressure loads.
The system stiffness matrix will have a degeneracy equal in

number to the rigid-body modes of the structure. On application
of adequate boundary conditions, rigid-body modes and the
degeneracy of the stiffness matrix are removed and hence, the
solution of the system of equations becomes possible (see
Application of boundary conditions on page 52 for a discussion
on the application of boundary conditions). The system of
equilibrium equations after the application of the boundary
conditions becomes
ˆ ˆ
[ K ] s { X }s = { R }s (8)
and a unique solution { X } s may be obtained. Once the nodal

displacements are obtained, the element stresses and nodal
forces may be computed. For the elements at the boundaries
these nodal forces will be in equilibrium with the reaction forces.
Friday 28 April 2000 – 14:43

Implementation S-Frame uses a skyline implementation (see Skyline storage

details scheme on page 48 for a discussion on the skyline storage
scheme) of the Crout reduction (Bathe (1982), Gerald (1980),
Press et. al. (1990)) to factorize the system stiffness matrix to
ˆ ˆ T
[ K ]s = [ L ] [ D ] [ L ] (9)
The stiffness matrix factorization, a computationally intensive

process, is performed only once during a linear static analysis
regardless of the number of load cases and load combinations.
Therefore, for each load case and/or load combination i, S-Frame
needs only forward reduce and back substitute for the
corresponding system load vector, i.e.
ˆ T ˆ
[ L ] [ D ] [ L ] = { R }s i (10)
During matrix factorization, the skyline solver performs the

following checks:
1. If a pivot is less than zero (system matrix is negative definite),

then it issues the SL001 exception. In such a case, the Crout
factorization can not be performed, therefore, it will abort
the factorization.
2. If a pivot is identically zero (system matrix is positive semi-

definite), then it issues the SL003 exception, it replaces the
zero pivot with a very small pivot and continues the
factorization.
3. If a pivot is greater than zero but less than tstval, then it

issues the SL002 exception. In such a case, the system matrix
may be ill-conditioned and results may be meaningless. The
Crout factorization will continue. The user should ensure
that the results from the analysis are physically plausible. The
value of the test variable tstval is calculated as follows:
– 31
tst ( dv )al = Max ( 7.888 ×10 , ( ftr1 ftr2 )maxpvt ) (11)
where
maxpvt maximum pivot encountered.
– 13
ftr1 1 ×10 (default)
– 23
ftr2 1 ×10 (used if support settlements are present
in the model)
Free-vibration The purpose of the free-vibration analysis is to determine the

natural frequencies and corresponding mode shapes of vibration.
analysis
This information is needed for instance for seismic analysis. The
following assumptions are made:
1. Stiffness and inertial effects of the structure are time

independent and hence the free-vibration motion is simple
harmonic.
2. Damping effects are ignored.
3. No forces, displacements, pressures or temperature effects

are applied to the structure.
Analysis description The equilibrium equations for the free-vibration analysis are
··
[ M ] s { X } s + [ K ]s { X } s = { 0 } (12)
In a free-vibration analysis the motion is simple harmonic and

may be expressed as
{ X }s = { φ } sin ( ωt + θ ) (13)
The acceleration vector may be obtained by taking the second

time derivative of Eq. 13
·· 2
{ X } s = – ω { φ } sin ( ωt + θ ) (14)
Substituting Eqs. 13 and 14 into Eq. 12 yields the following

expression
2
– ω [ M ]s { φ } sin ( ωt + θ ) + [ K ] s { φ } sin ( ωt + θ ) = { 0 } (15)
Friday 28 April 2000 – 14:43

which may be re-written as
2
( [ K ]s – ω [ M ]s ) { φ } = { 0 } (16)
since the sine term is arbitrary and is not equal to zero at all
times. Using Cramer’s rule (Anton (1977)), the solution of the
above equation can be shown to be of the form
{0}
{ φ } = ------------------------------------
2
- (17)
[ K ] s – ω [ M ]s
Hence, a non-trivial solution is possible only when the

determinant vanishes, i.e.
2
[ K ]s – ω [ M ]s = 0 (18)
which is the basic statement for free vibration problems. The

above equation is the frequency equation of the system and its
solution, after applying the necessary boundary conditions (see
Application of boundary conditions on page 52 for a discussion
on the application of boundary conditions), corresponds to the
solution of the general eigenvalue problem
ˆ 2
( [ K ]s – ωi [ M ] s ) { φi } = { 0 } (19)
for the numdof eigenpairs
ωi , { φi } (20)
Implementation There are numerous methods for solving the general eigenvalue
details problem stated in Eq. 19. For small systems (a few hundred
degrees of freedom) some of the most popular methods include
the Jacobi, Given’s and Householder methods (Petyt (1990)). For
larger systems, methods of reducing the number of degrees of
freedom have been developed and are described by Bathe (1982)
and Petyt (1990). However, for very large systems
aforementioned methods are inefficient. The eigenvalue solver
used in S-Frame is the skyline implementation (see Skyline
storage scheme on page 48 for a discussion on the skyline
storage scheme) of the Subspace Iteration method (Bathe
(1982)).
Subspace iteration is a very effective method of determining the

numnatfrq lowest eigenvalues and corresponding eigenvectors of
Eq. 19 simultaneously. The subspace iteration procedure is as
follows:
1. The number of eigenpairs (numnatfrq) to be extracted is set

to the user-specified number of eigenvalues. Note that
numnatfrq will normally be much smaller than the number
of degrees of freedom (i.e. numnatfrq << numdof). If this
condition is not met, the solver may issue the SI004 and/or
JD001 exceptions. For the current implementation of the
subspace iteration, it is recommended that numnatfrq be
below 100. For the case when the model is using the
consistent mass formulation (i.e. a non-zero density has been
specified for all the elements in the structure), the user may
easily calculate the structure’s number of degrees of freedom
by the following equation:
numdof = 6 numnod – numbc (21)

However, in the case where lumped masses are used to model
the inertia of the structure, the number of degrees of
freedom is equal to the number of lumped masses times the
degrees of freedom per lumped mass present in the structure.
2. Select a starting matrix [ A ]1 of dimension (numdof x

numitrvec) where numitrvec > numnatfrq and is given by
numitrvec = Min ( 2 numnatfrq, numnatfrq + 8 ) (22)

The columns of matrix [ A ]1 are determined as follows: The
first column of the product [ M ] s [ A ] 1 consists of the diagonal
terms of the system mass matrix [ M ] s . The other columns in
[ M ] s [ A ] 1 are unit vectors each with unit entries at the position
of the smallest ratio Kii ⁄ Mii where Kii and Mii are the diagonal
terms of the system stiffness and mass matrices respectively.
3. Perform the following operations:
(i) For i = 1, 2,... solve for the reduced matrix [ Ar ] i + 1 from

equation
ˆ
[ K ]s [ A r ] i + 1 = [ M ] s [ A ] i (23)
Friday 28 April 2000 – 14:43

Above system of equations is solved using the skyline

solver described in Implementation details on page 14.
(ii) Calculate the reduced stiffness and mass matrices
T ˆ
[ Kr ]i + 1 = [ Ar ] i + 1 [ K ]s [ Ar ]i + 1 (24)
T
[ Mr ]i + 1 = [ Ar ] i + 1 [ M ] s [ Ar ] i + 1 (25)
(iii)Solve the reduced general eigenvalue problem
[ Kr ]i + 1 [ Φ ]i + 1 = [ Mr ] [ Φ] i + 1 [ Λ ] (26)
For numnatfrq eigenvalues [ Λ ] i + 1 and eigenvectors
[ Φ ] i + 1 . This reduced eigenvalue problem may be solved
by any of the many methods available for the solution of
the general eigenvalue problem. Noting that the reduced
stiffness and mass matrices ( [ K r ] , [ Mr ] ) tend toward a
diagonal form as the number of subspace iterations
increases, the generalized Jacobi method is the most
effective choice.
(iv)Calculate an improved approximation to the eigenvectors
of the original system using
[ A ]i + i = [ A r ]i + 1 [ Φ ]i + 1 (27)
The eigenvalues [ Λ ] i + 1 and eigenvectors [ Φ ] i + 1 converge
to the lowest eigenvalues and eigenvectors of Eq. 19 as
i → ∞ . Convergency is achieved and this process is
terminated when the user-specified tolerance is met for
all numnatfrq eigenvalues. It should be noted that the
tolerance with which the lower eigenvalues are extracted
is smaller than the tolerance with which the higher
eigenvalues are extracted. A user-specified tolerance of
0.001 is generally recommended. If convergence is not
achieved within the user-specified maximum number of
iterations, the solver will issue the SI001 exception. The
recommended user-specified number of iterations is
around 12.
4. Use the Sturm Sequence check (Jennings (1977)) to verify

that no eigenvalues have been missed.
Shifting Procedure Above procedure for extracting the natural frequencies and
mode shapes of vibration will fail if the system stiffness matrix
ˆ
[ K ]s is singular. This case arises with not-fully-supported
structures commonly found in space applications. In order to

allow the user to perform a free-vibration analysis of such
structures, S-Frame uses the shifting procedure (Bathe (1982)).
In the solution of Eq. 19, S-Frame performs a shift λ on [ Kˆ ] s by
calculating
)
[ K ]s = [ K ] s – λ [ M ]s (28)
and then it considers the eigenvalue problem

)
2
( [ K ]s – ϖ i [ M ]s ) { φi } = 0 (29)
which yields the same eigenvectors { φ i } as the original

eigenvalue problem of Eq. 19 with its eigenvalues related to those
of Eq. 19 by
2 2
ϖi = ωi –λ (30)
The shift factor entered by the user, is used to calculate the

amount of shifting using the following relation
ˆ
RMS ( [ K ] s )
λ = shift factor --------------------------- (31)
RMS ( [ M ] s )
where RMS([ ]) is the root-mean-square of the diagonal

components of the corresponding matrix. Typically, a shift factor
of magnitude around unity is recommended for problems with
rigid-body motion.
Closely packed When the user suspects that the eigenvalues of a system are
eigenvalues closely packed, then he/she should run several analyses using
progressively smaller tolerance and larger number of iterations
to ensure that no eigenvalues have been missed. In addition,
shifting may be used to accelerate the convergence of closely
packed eigenvalues.
Mode shape If the value of one of the elements of a mode shape vector { φ i } is
normalization assigned a specified value, then the remaining elements are
determined uniquely. The process of scaling a natural mode so
Friday 28 April 2000 – 14:43

that each of its elements has a unique value is called

normalization, and the resulting modal vectors are called normal
modes (Craig (1981)).
Denoting by { ψ i } a mode that has been scaled to make its

amplitude unique, and assuming { ψ i } to be “dimensionless”, that
is, an arbitrary modal vector corresponding to ωi can be written
in the form
{ φ }i = λ i { ψ }i (32)
where λ i is a scaling factor whose units are such that

T
{ ψ } i [ M ] s { ψ }i has the dimensions of mass.(a)
S-Frame scales mode shapes so that the generalized mass defined

by
T
m i = { ψ } i [ M ]s { ψ } i (33)
has the value of unity.(b)
The generalized stiffness for the ith mode is defined as
T
ki = { ψ }i [ M ] s { ψ }i (34)
In the case where the mode shapes are normalized with respect
to the mass matrix (Eq. 33), generalized stiffnesses are given by
2
ki = ωi (35)
(a) The reason for this particular definition of dimensionless is that vectors { ψ }i may
contain a mixture of types of coordinates, for example, translations and rotations.
Thus, it would not be possible to simultaneously make all components of { ψ }i
dimensionless in the usual sense of the word.
(b) As previously noted, it is convenient to scale { ψ }i , so that the product
T
{ ψ } i [ M ] s { ψ } i has the units of mass, therefore, m i = 1Kg (or 1 Slug).
In addition to the natural frequencies and mode shapes of

vibration, S-Frame computes for each mode i the modal
j
participation factors( jµi ) and modal masses ( mˆ i ) in the three
global directions (j = X, Y, Z) defined, respectively, as
T j
j – { ψ } [ M ] { I }-
µi = ----------------------------------
T
(36)
{ψ} [M]{ψ}
T j 2
jˆ ({ψ } [M] {I })
mi = --------------------------------------
T
- (37)
{ ψ} [M]{ ψ}
where {j I } is a unit vector in the jth global direction which is

expanded to the system dimension. For instance Y{ I } is given by
010000 0 10000 010 000 … (38)
Since S-Frame scales the mode shapes with respect to the mass
matrix, the expressions for the modal participation factors and
modal masses simplify to
j T j
µi = – { ψ } [ M ] { I } (39)
jˆ T j 2
mi = ( { ψ } [ M ] { I } ) (40)
Using the modal masses, S-Frame also computes the percentage

of the mass in each of the three global directions (j = X, Y, Z)
using the following expression
numnatfrq
jˆ
∑ mi
j i=1
% M = -----------------------------
j
- (41)
[ M] {I }
Stressed free- The purpose of the stressed free-vibration analysis is to

determine the natural frequencies and corresponding mode
vibration analysis
shapes of vibration while taking into consideration membrane
effects due to time-independent loads under the following
assumptions:
1. Stiffness and inertial effects of the structure are constant.
Friday 28 April 2000 – 14:43

2. Damping effects are ignored.
3. No time-dependent forces, displacements, pressures, or

temperature effects are applied to the structure.
4. Time-independent forces, displacements, pressures, or

temperature effects may be applied to the structure.
Analysis description The stressed free-vibration analysis requires a two phase analysis
procedure. In the first phase, the system equilibrium equations
for the linear static analysis (Eq. 7) are used to solve for the { X } s
for the user-selected load case or load combination. Once the
nodal displacements are obtained, the element membrane forces
are calculated and used to form the element geometric stiffness
matrices [ Kg ]e . In general, compressive membrane forces tend
to reduce the stiffness of an element, while tensile membrane
forces cause a corresponding increase of stiffness. In the second
phase of the analysis, the equilibrium equations for the stressed
free-vibration analysis are considered, namely
··
[ M ]s{ X }s + [ [ K ]s + [ Kg ]s ] { X } s = { 0 } (42)
where the system geometric stiffness matrix [ Kg ] s is obtained by

assembly of the element geometric stiffness matrices [ Kg ] e using
standard matrix procedures. The geometric stiffness matrix
represents the effect of membrane forces on the stiffness of the
elements. Hence, by including it in the equilibrium equations 42,
we are in essence including the effects of the membrane forces on
the stiffness of the structure. Following the same procedure as in
Analysis description on page 15, it can be shown that the
frequency equation of the stressed system is
2
[ [ K ] s + [ Kg ] ] – ωi [ M ]s = 0 (43)
and its solution, after applying the necessary boundary

conditions (see Application of boundary conditions on page 52
for a discussion on the application of boundary conditions),
corresponds to the solution of the general eigenvalue problem
ˆ 2
( [ [ K ] s + [ K g ]s ] – ω i [ M ] s ) { φ i } = 0 (44)
for the numdof eigenpairs
ω i , { φi } (45)
In general, the eigenvalues ωi (hence the natural frequencies of

oscillation) decrease in the presence of compressive membrane
forces, whereas they increase in the presence of tensile
membrane forces. In the case where the structure is critically
loaded, the lowest eigenvalue becomes zero. This is due to the
singularity of the stiffness matrix induced by the geometric
stiffness matrix. This is in fact one of the methods used to predict
buckling loads by frequency measurements for varying
membrane loads. Further details are provided on the methods
for determining buckling loads in Linear buckling analysis on
page 31.
Implementation The implementation of the stressed free-vibration analysis is a

details combination of the descriptions in Implementation details on
page 14 and Implementation details on page 16. In the current
version of S-Frame, the geometric stiffness matrices [ K g ] e of
beams and shells are formed.
Linear dynamic The purpose of the linear dynamic analysis is to determine the
response of a structure to arbitrarily time-varying loads under
analysis
the following assumptions:
1. Stiffness and inertial effects of the structure do not depend

on time.
2. Initial conditions are known.
3. No gyroscopic or Coriolis effects are included
Analysis description The system equilibrium equations for the linear dynamic
analysis are
·· ·
[ M ] s { X } s + [ C ] s { X } s + [ K ] s { X } s = { R }s (46)
The above system of equations is a system of linear second order

differential equations with constant coefficients. Three categories
of methods exist for the solution of such system of equations:
Friday 28 April 2000 – 14:43

direct explicit, direct implicit, and indirect. Direct methods do

not require any transformation of the system equilibrium
equations prior to solution. Indirect methods require the
transformation of the system equilibrium equations to a more
convenient form prior to solution (Modal Analysis). Direct
methods enforce equilibrium which includes the external forces,
·· ·
inertia forces ( [ M ] s { X } s ) and damping forces ( [ C ]s { X }s ) at discreet
time intervals. Methods that consider equilibrium at time t are
classified as explicit (Central Difference), whereas methods that
consider equilibrium at time t + ∆t are classified as implicit
(Newmark, Wilson, Linear Acceleration). S-Frame uses the
Newmark time integration method for the solution of Eq. 46.
The cost of dynamic analysis is proportional to the number of
steps in the analysis. Therefore, the choice of the step size
becomes very important. This point will be discussed further in
the following sections. Dropping the subscript s from all system
vectors for brevity, the difference equations suggested by
Newmark may be expressed as
· ·· ··
{ X } t + ∆t = { X }t + ∆t ( ( 1 – δ ) { X } t + δ { X }t + ∆t ) (47)
and
{ X } t + ∆t = { X } t + ∆t { X } t + ∆t 2  -- – α { X }t + α { X }t + ∆t

· 1 ·· ··
(48)
2
If the values of 1/2 and 1/6 are used for δ and α , respectively,
then above equations reduce to the equations used in the linear-
acceleration method. First we consider the system of second-
order differential equations at time t = t+∆t
·· ·
[ M ] s { X } t + ∆t + [ C ]s { X } t + ∆t + [ K ] s { X }t + ∆t = { R } t + ∆t (49)
from which we intent to solve for the displacement { X } t + ∆t . To do

·· ·
so, we must first obtain expressions for { X } t + ∆t and { X } t + ∆t in
terms of quantities at time t and { X } t + ∆t only. From Eq. 48 we
obtain
{ X } t + ∆t = -----------2- ( { X } t + ∆t – { X }t ) – ---------- { X }t –  ------- – 1 { X } t

·· 1 1 · 1 ··
 2α 
(50)
α∆t α∆t
which upon substitution in Eq. 47 yields
· · ·· ··
{ X }t + ∆t = { X } t + ∆t ( ( 1 – δ ) { X }t + δ { X } t + ∆t ) (51)
Substituting Eqs. 50 and 51 into the system of equations 49 and

rearranging, we obtain
[ K eq ]s { Y } = { R } t + ∆t + { R eq } t (52)
t + ∆t
where
1 δ
[ K eq ]s = -----------2- [ M ] s + ---------- [ C ] s + [ K ] s (53)
α∆t α∆t
and
{ R eq } t = [ M ] s  -----------2- { X }t + ---------- { X } t +  ------- – 1 { X } t +

1 1 · 1 ··
 α∆t α∆t  2α  
(54)
δ δ δ
[ C ] s  ---------- { X }t +  --- – 1 { X } t + ∆t  ------- – 1 { X } t
· ··
α∆t α 2α
The system of equations in 52 can then be solved to obtain the

displacements at t = t+∆t. Having done that, Eqs. 50 and 51 may
be used to obtain the velocities and accelerations at time t = t+∆t.
Since equilibrium is considered at t = t+∆t, the Newmark
method is an implicit time-integration method.
Implementation This section describes the Newmark algorithm as implemented

details of Newmark in S-Frame.
algorithm
Preliminary 1. Form the system matrices and apply appropriate boundary

calculations conditions to obtain [ Mˆ ] s and [ Kˆ ] s .
2. Select values for the Rayleigh damping parameters α and β

(see Damping effects on page 57 for more details), and form
the system damping matrix
ˆ ˆ ˆ
[ C ] s = α [ M ] s + β [ K ]s (55)
Damping effects on page 57 gives details on how to select the
damping parameters α and β .
Friday 28 April 2000 – 14:43

3. Select values for the Newmark parameters α and δ. If no

physical damping is present in the structure (i.e. Rayleigh
damping parameters α and β are zero) the values 0.2525 and
0.5050 should be used respectively. Otherwise, the values
0.25 and 0.50 should be used instead.
4. Select a small enough time-step ∆t. In general the value

0.00001 s should be used for the initial time-step size of the
variable time-step algorithm. In the case of the constant
time-step algorithm, the time-step size should be selected
using the guidelines given in Selection of constant time-step
size on page 28.
5. For the variable time-step algorithm, select a value for the Jcr
parameter. A description of this parameter is given in
Implementation details of adaptive time-step size
algorithm on page 27.
6. For the variable time-step algorithm, select values for the

minimum and maximum time-step sizes using the guidelines
of Selection of minimum and maximum time-step sizes on
page 29.
7. Using the initial displacement and velocity vectors,

determine the initial accelerations using equilibrium at time
to
ˆ ·· ˆ ˆ ·
[ M ] s { X }to = { R } t – [ K ] s { X } t – [ C ] s { X } to (56)
o o
For each time-step 1. Form [ K eq ]s
1 ˆ δ ˆ ˆ
[ K eq ]s = -----------2- [ M ] s + ---------- [ C ]s + [ K ] s (57)
α∆t α∆t
2. Factorize [ Keq ]s
T
[ K eq ]s = [ L ] [ D ] [ L ] (58)
3. Form { Req } t
{ R eq } t = [ M ] s  -----------2- { X }t + ---------- { X } t +  ------- – 1 { X t } +

ˆ 1 1 · 1 ··
 α∆t α∆t  2α  
(59)
δ δ δ
[ C ] s  ---------- { X }t +  --- – 1 { X } t + ∆t  ------- – 1 { X } t
ˆ · ··
 α∆t α   2α  
4. Back-substitute and forward reduce for the displacements at

t = t+∆t
T ˆ
[ L ] [ D ] [ L ] { X }t + ∆t = { R } t + ∆t + { R eq } t (60)
5. Solve for the accelerations and velocities at t = t+∆t using

Eqs. 50 and 51 respectively.
6. For the variable time-step size algorithm, adapt the time-step

size according to the procedure described in
Implementation details of adaptive time-step size
algorithm on page 27.
Implementation The criterion used for the selection of a time-step size is the
details of adaptive magnitude of the rate of change of acceleration, known as jerk.
time-step size Jerk is approximated by
algorithm ·· ··
{ X } n – { X }n – 1
{ J } n = ------------------------------------- (61)
∆t n
where the most recent results available are via the nth integration
step using a time-step size of ∆tn. Next we determine the
maximum component of the jerk vector, which we denote by Jˆ n .
The maximum jerk is then normalized as follows
ˆ
Jn = J n ∆t n3 (62)
The user of this procedure selects a value for the critical

normalized jerk, denoted by Jcr, which corresponds to the largest
value of normalized jerk allowed at the n+1 integration step, that
is
ˆ
J cr = Jn + 1 = J n + 1 ∆t n3 + 1 (63)
Friday 28 April 2000 – 14:43

For small time-step sizes, the above expression may be

approximated by
ˆ
J cr = Jn ∆t n3+ 1 (64)
Now substituting for ˆJn from Eq. 62 and solving for ∆tn+1, Eq. 64
yields an expression for the estimation of the next time-step size,
namely
Jcr 1 / 3
∆t n + 1 = ∆t n  ------ (65)
Jn
For structural dynamics problems, the user-specified value of Jcr

should be in the neighborhood of 0.001. From a number of tests
we have performed, we found that a 30% reduction in the
number of time steps is normal in structural dynamics when
using the variable time-step size algorithm. The adaptive time-
step algorithm performs numerous checks before selecting an
appropriate time-step size. One of these checks is to ensure that
the time-step size does not increase or decrease too rapidly (a
known cause of numerically induced transients). When this
occurs, the DE005 exception is issued. If this exception occurs at
the beginning of the analysis only, then the value of the user-
defined initial time-step size should be adjusted accordingly. If it
occurs at the end of the analysis it may be ignored. If, however, it
occurs throughout the analysis, a smaller value for the user-
specified value of Jcr is recommended.
Selection of constant It was noted (Tabarrok and Stylianou (1985)) that the response in
time-step size modes with a small ratio Ω ⁄ ωi is essentially static, and the
response in modes with a large ratio Ω ⁄ ωi is negligible (where Ω
is the highest frequency component of the forcing function).
Therefore, realizing that the static response is directly included in
the time integration of the system equilibrium equations 46, we
need only consider response contributions from modes for
which the ratio Ω ⁄ ωi is less or equal to 1 ⁄ 4 . Two important
practical considerations follow from the above discussion.
Firstly, the time integration need only be accurately performed
for modes with frequencies up to 4Ω , which will be denoted by
*
Ω from here on. And secondly, the finite-element model used
need only accurately represent modes with frequencies up to Ω * .
As a rule of thumb we recommend that the user-selected
constant time step size satisfies the following condition
2π
∆t cr ≤ ------------* (66)
20Ω
Selection of As a rule of thumb, we recommend that the user-selected

minimum and minimum and maximum time-step sizes satisfy the following
maximum time-step conditions, respectively
sizes
∆t cr
--------- ≥ ∆tmin ≥ 1.0 ×10– 8 (67)
10.0
∆t max ≥ 100.0∆t cr (68)
where ∆t cr is estimated using Eq. 66.
P-delta analysis The purpose of the P-delta analysis is to determine the

1. The effect of axial loads (membrane forces) on the stiffness

of the structure is considered.
2. Stiffness effects and applied loads do not depend on time.

(i.e. dead loads).

Analysis description The bending moments of beams can be of two types:
1. Primary bending moments – caused by moments applied at

the ends of the beam or by transverse loading along the span
of the beam.
Friday 28 April 2000 – 14:43

2. Secondary bending moments – caused by the axial force

acting through the lateral displacement of the beam.
The so-called P-δ moments are moments caused by the axial

force acting through the lateral displacement of the beam relative
to its chord (Fig. 1.a), whereas the P-∆ moments are moments
caused by the axial force acting through the relative lateral
displacement of the two ends of the beam (Fig. 1.b). The effect of
axial force on the primary bending moments of slender beams
and shells may be significant. The procedure in which these
secondary effects are accounted for in a static analysis in referred
to as a P-delta analysis.
δ
P P
(a)
∆
(b)
Figure 1: (a) P-δ effect, (b) P-∆ effect.

The P-delta analysis, specifically the Two-Cycle Iterative Method
(Chen and Lui (1991)) requires a two pass analysis procedure. In
the first phase, the system equilibrium equations for the linear
static analysis (Eq. 7) are solved for the nodal displacements { X } s
for all of the user-selected load cases and load combinations.
Once the nodal displacements are obtained, the element
membrane forces are calculated and used to form the element
geometric stiffness matrices [ Kg ] e . The beam geometric stiffness
matrix includes both the P-δ and P-∆ effects. As mentioned

earlier, compressive membrane forces tend to reduce the stiffness
of an element, while tensile membrane forces cause a
corresponding increase of stiffness. In the second phase of the
analysis, the equilibrium equations of the P-delta analysis are
considered, namely
( [ K ]s + [ Kg ]s ) { X }s = { R }s (69)
where the system geometric stiffness matrix [ K g ] s is obtained by

assembly of the element geometric stiffness matrices [ Kg ]e .
Accordingly, Eqs. 69 are solved for the nodal displacements { X } s
for all of the user-selected load cases and load combinations.
Once nodal displacements are obtained, the element stresses and
nodal forces may be computed. For the elements at the
boundaries, these nodal forces will be in equilibrium with the
reactions.
Implementation S-Frame’s implementation of the two-cycle iterative P-delta

details analysis is an approximation to the full non-linear analysis. For
the solution of the system of equations, the skyline solver as
described in Implementation details on page 14 is used. In the
current version of S-Frame, the geometric stiffness matrices [ Kg ] e
of beams and shells are formed.
Linear buckling The purpose of the linear buckling analysis is to determine the
critical load factors and corresponding buckling mode shapes
analysis
under the following assumptions:
1. Stiffness of the structure does not change with time.
2. The structure’s response to load is unbounded (i.e. with a

horizontal force-deflection curve).
Analysis description The linear buckling analysis requires a two phase analysis
procedure. In the first phase, the system equilibrium equations
for the linear static analysis (Eq. 7) are used to solve for the
displacements { X } s for the user-selected load case or load
combination. Once the nodal displacements are obtained, the
element membrane forces are calculated and used to form the
Friday 28 April 2000 – 14:43

element geometric stiffness matrices [ Kg ] e . As mentioned earlier,

compressive membrane forces tend to reduce the stiffness of an
element, while tensile membrane forces cause a corresponding
increase of stiffness. In the second phase of the analysis, the
equilibrium equations for the linear buckling analysis are
considered, namely
[ [ K ]s + [ Kg ]s ] { X }s = { R }s (70)
where the system geometric stiffness matrix [ Kg ] s is obtained by

assembly of the element geometric stiffness matrices [ Kg ] e , using
standard matrix procedures. The geometric stiffness matrix
represents the effect of membrane forces on the stiffness of the
elements. Hence, by including it in the equilibrium equations 70,
we are in essence including the effects of the membrane forces on
the stiffness of the structure. In order to evaluate the critical
buckling loading of the structure let us express the geometric
stiffness matrix in terms of a reference geometric stiffness matrix
multiplied by a load factor λ . Thus
[ Kg ]s = λ [ K g ] (71)
o s
where [ Kg o ]s is the geometric stiffness matrix of the structure

under the reference loading (in general { R } s is used as the
reference loading). Therefore, the loading of the structure is
proportional to λ and its relative distribution is constant (Cook
et. al. (1989)). Substituting Eq. 71 in Eq. 70 the equilibrium
equations become
( [ K ] s + λ [ K g ] ) { X }s = { R } s (72)
o s
For instability to occur, the equivalent stiffness matrix
[ K eq ]s = [ K ] s + λ [ K g ] (73)
o s
must become singular. Hence, the most obvious way for solving
the above problem, is to increase the load factor, step by step, and
at each load level check the singularity of matrix [ Keq ] s by
evaluating its determinant. Three possibilities exist:
[ K eq ]s < 0 if ( λ > λ cr )
[ K eq ] s = 0 if ( λ = λcr ) (74)
[ K eq ] s > 0 if ( λ < λ cr )
where λ cr is the critical load factor that if applied loads ( { R } s )

are multiplied by, will yield a loading configuration that will
render the structure unstable (i.e. it will cause the onset of
buckling). The first one corresponds to unstable equilibrium, the
second possibility corresponds to the case where the structure is
on the verge of instability, and the third one corresponds to a
stable structure. The value of λ is then a measure of the closeness
of the actual loading to the loading necessary for the structure to
buckle. This problem can also be viewed as an eigenvalue
problem (for more details see Cook (1989), Tabarrok and
Stylianou (1985))
2
[ K ]s + λ i [ Kg ] = 0 (75)
o s
which can be identified as the general eigenvalue problem

described in Analysis description on page 15 for the free-
vibration analysis. Now all that is necessary, is to determine the
smallest eigenvalue and corresponding eigenvector. Generally,
there is interest only in the first eigenpair
λ cr , { φ 1 } (76)
where, as mentioned earlier, λ cr is the critical load factor and

{ φ1 } the corresponding buckling mode shape.
Friday 28 April 2000 – 14:43

It should be noted that the buckling of the structure as treated

here, does not require or imply buckling of each element. That is,
at the structure’s critical loading, the following is true for an
element
2
( [ K ] e + λ cr [ K g ] ≥ 0 ) (77)
o e
where equality prevails in the highest buckling mode or for a

one-element structure.
Implementation The implementation of the linear buckling analysis is a

details combination of the descriptions in Implementation details on
page 14 and Implementation details on page 16. In the current
version of S-Frame, the geometric stiffness matrices [ Kg ] e of
beams and shells are formed. In addition, it is recommended that
only the first few buckling modes be computed in a linear
buckling analysis. Specifically, for small structures no more than
the first couple mode shapes should be requested.
Buckling - Stressed It is instructive to consider the three possibilities described by

Free Vibration - P- Eqs. 74 in the context of two other analysis types, namely
delta Stressed Free Vibration and P-delta analyses and compare to the
Buckling analysis.
Buckling Analysis [ K eq ] s < 0 if ( λ > λ cr ) : the absolute value of λ cr will be less than 1
indicating that the structure has already buckled.
[ K eq ] s = 0 if ( λ = λ cr ) : the absolute value of λ cr will be identically

1 indicating that the structure is on the verge of instability.
[ K eq ] s > 0 if ( λ < λ cr ) the absolute value of λ cr will be greater

than 1 indicating that the structure is stable
Stressed Free Vibration [ K eq ] s < 0 if ( λ > λ cr ) : the analysis will fail since at least one of
Analysis
the eigenvalues is negative.
[ K eq ] s = 0 if ( λ = λ cr ) the analysis may be successful with a zero

lowest natural frequency.
[ K eq ] s > 0 if ( λ < λ cr ) the analysis will be successful with all

positive eigenvalues.
P-delta Analysis [ K eq ] s < 0 if ( λ > λ cr ) : the analysis will fail since at least one of
the pivots is negative.
[ K eq ] s = 0 if ( λ = λ cr ) at least one of the pivots will be very close

to zero. If it is zero or positive then the analysis will not fail, but
the displacements will be physically large. If it is negative the
analysis will fail.
[ K eq ] s > 0 if ( λ < λ cr ) the analysis will be successful with all the

pivots greater than zero.
Response The purpose of the response spectrum analysis is to determine

the stresses due to a response spectrum under the following
spectrum analysis
assumptions:
1. The structure is linear.
2. Structure is excited at all support nodes by the same random

spectrum.
3. Spectrum’s direction and frequency content are known.
Friday 28 April 2000 – 14:43

Analysis description In the case of transient excitation at the supports of a structure, a

detailed time history of the response is not always required. The
response spectrum analysis estimates the peak response, such as
maximum displacement, without the need for computing the
time history of the response.
Y’
{W}s
Y
{X}s X’
{Z}s
Z X
Figure 2: Ground motion.

The system equilibrium equations associated with the response
spectrum analysis are
·· · ··
[ M ] s { W } s + [ C ]s { W } s + [ K ] s { W } s = –[ M ]s { Z ( t ) } s (78)
where the relative displacement vector { W } s and ground

displacement vector { Z } s (a given function of time t) are shown
in Fig. 2 together with the absolute displacement vector { X } s .
The steps required for a response spectrum analysis are as

follows:
1. Perform a free-vibration analysis (ignoring damping) to

compute the first numnatfrq eigenpairs ωi , { φ i }.
2. Form the normalized modal matrix (see Eq. 33)
[ Ψ ] = [ { ψ 1 }, { ψ 2 }, { ψ 3 }, …, { ψ numnatfrq } ] (79)
3. Using the normalized mode-shape matrix transform the

system equilibrium equations 78
T ·· ··
[ Ψ ] [ M ]s [ Ψ ] { η }s = [ m ]s{ η }s
T · ·
[ Ψ ] [ C ]s [ Ψ ] { η } s = [ c ] s { η } s
(80)
T
[ Ψ ] [ K ]s[ Ψ ] { η }s = [ k ]s { η }s
T ·· T j ··
–[ Ψ ] [ M ]s { Z ( t ) }s = –[ Ψ ] [ M ]s { I } z ( t )
where we have separated the base acceleration vector into a

constant vector which indicates the global direction of the
base motion ( {j I } , where j = X, Y, Z), and a scalar function
which gives the ground acceleration history in time ( z··( t ) ). If
the mode shapes are normalized according to Eq. 33, then
the system equilibrium equations will be decoupled so that
for the ith generalized degree of freedom the equilibrium
equation becomes
·· · 2 j ··
η i + 2 ξ i ω i η i + ω i η i = µi z ( t ) (81)
where ξi is the modal damping ratio defined in Damping
effects on page 57, and jµi is the participation factor of mode
i in direction j and is defined as
j T j
µi = – { ψ } i [ M ] s { I } (82)
Friday 28 April 2000 – 14:43

4. The maximum relative displacement contribution of mode i

due to the ground motion in the jth direction can be
expressed in terms of response spectra as
j
j max µi
{ w }i = { ψ }i ----- S v ( ω i, ξ i ) (83)
ωi
where S v ( ω i, ξ i ) is the spectral pseudo-velocity (Clough and

Penzien (1975), Craig (1981)).
5. Combine the maximum relative displacement contribution

of mode i due to the ground motion in all directions using
one of the following methods:
Absolute Sum of Modal Maxima

numdir
max j max
wi = ∑ wi (84)
j=1
Square Root of the Sum of the Squares

numdir
max j max 2
{ w }i = ∑ ( { w }i ) (85)
j=1
6. Combine the maximum modal contributions using one of

the following methods:
Absolute Sum of Modal Maxima

numnatfrq
max max
W = ∑ wi (86)
i=1
Square Root of the Sum of the Squares

numnatfrq
max max 2
{W} = ∑ ( { w }i ) (87)
i=1
Complete Quadratic Combination

numnatfrq
max max max
{ W} = ∑ ∑ { w }i κ ik { w }k (88)
i=1 k=1
where
3⁄2
8 ξ i ξ k ( ξ i + rξ s )r
κ ik = ------------------------------------------------------------------------------------------------
- (89)
2 2 2 2 2 2
( 1 – r ) + 4ξ i ξ k r ( 1 + r ) + 4 ( ξ i + ξ k )r
ωi
r = ------ . (90)
ωk
Double Sum Combination Method

numnatfrq
max max max
{W} = ∑ ∑ κ ik { w }i { w }k (91)
i=1 k=1
where
1
κ ik = ----------------------------------------------------2 (92)
ω′ i – ω′ k
1 + -----------------------------------
ξ′i ω′ i – ξ′k ω′ k
2
ω′ i = ω i 1 – ξi (93)
2
ξ′ i = ξi + ---------- (94)
τe ωi
where τ e is the time duration of the white noise segment of

the earthquake excitation (represented by
extremely irregular accelerations of roughly equal
intensity).
Grouping Method
numnatfrq
max max max
{W} = ∑ ∑ κ ik { w }i { w }k (95)
i=1 k=1
where
ωk – ω i
1 if ----------------- ≤ 0.1
ωi
κ ik = (96)
ωk – ω i
0 if ----------------- > 0.1
ωi
Friday 28 April 2000 – 14:43

Implementation The implementation of the response spectrum analysis begins

details with a free-vibration analysis (Implementation details on page
16) and then a post-processing phase which uses equations 84 to
96 to estimate the maximum displacements, member forces,
reactions and base shears.
For the Complete Quadratic and Double Sum combination

methods, the current version of S-Frame uses a constant value of
modal damping ratio ξ i for all modes. This value is specified by
the user when choosing the modal combination method.
In a response spectrum analysis, it is important that the analyst

includes enough number of modes. If physically possible, at least
90% of the total mass of the structure should be recovered in the
directions of interest. Therefore, it is recommended that the
analyst perform several free-vibration analyses to determine how
many modes are necessary and then perform a response
spectrum analysis.
The current version of S-Frame performs the spatial

combination first and then the modal combination. The
implementation is given below:
Spatial combination for each mode
for each active direction
combine modal displacement contributions using

selected spatial combination method to get the total
modal displacement contributions
for each mode
for each element
compute the modal element force contributions

using the modal displacement contributions (not the
total)
combine the modal element force contributions

using selected spatial combination method to get the
total modal element force contributions
for each mode
compute the modal reaction contributions using the

modal displacement contributions (not the total)
combine the modal reaction contributions using selected

spatial combination method to get the total modal
reaction contributions
If more than one direction is active, the total modal

displacements, the element forces and the reaction contributions
will have been computed using either a power (square) or an
absolute value operator both of which are non-commutative
operators.
Modal combination for each mode
combine the total modal displacement contributions using

selected modal combination method to get the maximum
displacements
for each mode
combine the total modal element force contributions using

element forces
for each mode
combine the total modal reaction contributions using

reactions
If more than one mode is requested, then the maximum

displacements, the element forces and the reactions will have
been computed using one of the non-commutative RSA
combination methods. Therefore, due to the non-commutative
Friday 28 April 2000 – 14:43

nature of the spatial and modal combination methods, the

computed maximum quantities will differ if the modal
combination is performed before the spatial.
Loss of sign Because the maximum quantities (displacements, element forces

and reactions) are computed using statistical methods, all sign
information is lost. To remedy this situation, S-Frame assigns
signs to the maximum quantities as obtained from the dominant
modal contribution of the corresponding quantity in the
direction of the applied spectral curve.
Scaling S-Frame performs an additional scaling of all the computed

maximum quantities. The type of scaling factor is user specified
and it depends on whether the analysis is RSA or static.
If the analysis type is RSA, S-Frame scales all the maximum

quantities (displacements, reactions, element forces, base shears)
using the following scaling factors:
Tx Ty Tz
S x = ----- S y = ----- S z = ----- (97)
Rx Ry Rz
If the analysis type is Static, S-Frame scales only the dominant

mode quantities (displacements, reactions, element forces, base
shears), which then may be used to combine with static load
cases. The scaling factors used depend on the user selection and
are given by:
Tx Ty Tz
S x = ----- S y = ----- S z = ---- (98)
Fx Fy Fz
or
Rx Ry Rz
S x = ----- S y = ----- S z = ----- (99)
Fx Fy Fz
where
S i Scaling factor in direction i.
Ti Code base shear in direction i (user specified).
R i RSA Base Shear in direction i.
F i Dominant modal base shear in direction i.
Nonlinear static The purpose of the nonlinear static analysis is to determine the
analysis
1. Stiffness effects and applied loads may depend on

deformation.

(i.e. dead loads).

Analysis description The system equilibrium equations for the nonlinear static
analysis are
{ F }s – { R } s = { 0 } (100)
where the system load vector { R } s includes the contributions of

all externally applied nodal forces
1. Applied nodal loads
2. Loading due to static acceleration fields (such as gravity)
3. Element thermal/pressure loads
and vector { F }s includes the nodal forces that correspond to the

element stresses in this configuration. This vector is known as
the internal load vector.
Implementation The solution of the nonlinear Eq. 100 is carried out in an

details incremental step-by-step analysis with the total applied loads
divided into a number of load steps (n). Let us assume that the
Friday 28 April 2000 – 14:43

solution of Eq. 100 is known at load increment i (i.e. the system

displacement vector i{ X }s is known), and that the solution is
sought for load increment i+1, i.e.
i+1 i+1
{ F }s – { R }s = { 0 } (101)
Since the solution at load increment i is known, the internal load

vector at load increment i+1 may be written as
i+1 i
{ F }s = { F s } + ∆ { F s } (102)
where ∆ { Fs } is the change in the internal load vector

corresponding to the increment in displacements ∆ { Xs } from
load increment i to i+1. In this step-by-step analysis, the
incremental internal load vector is approximated by
i
∆ { Fs } = [ Kt ] ∆ { X s } (103)
where [ Kt ] is known as the tangent stiffness matrix. Substituting

Eqs. 102 and 103 into Eq. 101 we get the incremental equilibrium
equations
i i+1 i
[ Kt ] ∆ { Xs } = { R s } – { Fs } (104)
After solving for the incremental displacement vector, the

displacements at load increment i+1 can be estimated by
i+1 i
{ Xs } = { Xs } + ∆ { Xs } (105)
and one may then proceed to the next load increment. In

practice, however, this procedure will invariably be subject to
significant errors due to the assumption of Eq. 103, hence one
iterates until the solution of Eq. 101 is obtained with an
acceptable accuracy. Many methods are available for performing
these iterations. The next section describes one of the most
popular iteration methods used for the solution of nonlinear
equations, the Newton-Raphson method.
Newton-Raphson The most popular method for the solution of nonlinear

Method equations is the Newton-Raphson method (which is also used by
the current version of S-Frame). The Newton-Raphson iterations
(j) are performed at each load increment (i) using the following
equations
i j i+1 i+i j–1

[ Kt ] ∆ { X s } = { Rs } – { Fs }
(106)
i+1 j i+i j–1 j
{ Xs } = { Xs } + ∆ { Xs }
and at the beginning of the iterations we use the initial values

given by
i+1 0 i
{ Xs } = { Xs }
(107)
i+1 0 i
{ Fs } = { F }
The right-hand side of the first of Eqs. 106 may be interpreted as

an out-of-balance load vector. A load vector that is not yet
balanced by the internal load vector (element stresses). In
essence, the goal of the Newton-Raphson iterations is to
minimize this out-of-balance load vector.
Convergence Criteria Various convergence criteria may be used to determine when to

stop the Newton-Raphson iterations (j) in order to proceed to
the next load increment (i). Some of these criteria are based on
displacements, some on residual forces, and some on a
combination of both (using some energy measure). For a
detailed description of these criteria, please refer to Bathe (1982).
In this section we will describe the convergence criteria used by
S-Frame.
1. Absolute Displacement: The absolute value of each

displacement component of iteration vector ∆ { Xs } j is
compared to a user-specified tolerance. If it is less than the
user-specified tolerance, the Newton-Raphson iterations stop
and the solution proceeds to the next load increment (i).
2. Absolute Displacement and Rotation: The absolute value of

each displacement and rotation component of the iteration
vector ∆ { Xs } j is compared to a user-specified tolerance. If it is
Friday 28 April 2000 – 14:43

less than the user-specified tolerance, the Newton-Raphson

iterations stop and the solution proceeds to the next load
increment (i).
3. Relative Displacement: The ratio of the absolute value of

each displacement component of the iteration vector ∆ { Xs } j
to the increment vector ∆ { X s } is compared to a user-specified
tolerance ratio. If it is less than the user-specified tolerance,
the Newton-Raphson iterations stop and the solution
proceeds to the next load increment (i).
4. Relative Displacement and Rotation: The ratio of each

displacement and rotation component of the iteration vector
j
∆ { X s } to the increment vector ∆ { X s } is compared to a user-
specified tolerance ratio. If it is less than the user-specified
tolerance, the Newton-Raphson iterations stop and the
solution proceeds to the next load increment (i).
Friday 28 April 2000 – 14:43

S-Frame Theory Manual Chapter 4 : Other procedures : 48
4 Other procedures
This section describes the storage scheme as well as various
procedures used in S-Frame.
Skyline storage The S-Frame finite-element analysis engine uses the skyline
storage scheme for storing the system stiffness, damping and
scheme
mass matrices (Bathe (1982)) for all the supported types of
analyses. To illustrate this storage scheme, consider the
hypothetical stiffness matrix shown below (since this matrix is
assumed to be symmetric, only the upper half is shown).
k11 k 12 0 k 14 0 0 0 0
k 22 k 23 0 0 k 26 0 0
k 33 k 34 0 0 0 k 38
k 44 k 45 k 46 0 0
[K] = (108)
k 55 k 56 0 0
k 66 k67 0
k77 k 78
k 88
Therefore, to store a square symmetric matrix using a symmetric

storage scheme, where only the upper part of the matrix
(including the diagonal elements) is stored, requires storage for
n (n + 1 )
-------------------- = 36 (109)
2
49 : Chapter 4 : Other procedures S-Frame Theory Manual
elements, where n is the order of the square matrix. On the other

hand, if a symmetric banded storage scheme is used, the matrix
of Eq. 108 is stored as
k 11 k 12 0 k 14 0 0
k 22 k 23 0 0 k 26 0
k 33 k 34 0 0 0 k 38
k 44 k 45 k 46 0 0
[ K] = n
k 55 k 56 0 0 (110)
k 66 k 67 0
k 77 k 78
k 88
wh
hence, only
wh ( wh – 1 )
nw h – --------------------------- = 33 (111)
2
elements need to be stored, where wh is the half bandwidth. One

can easily see that the storage savings by using a symmetric
banded storage scheme can be significant if wh is much smaller
than n.
In the case of the skyline storage scheme, only the elements

between the diagonal and the last non-zero element of each
column are stored (including the diagonal). Therefore, the
matrix of Eq. 108 is stored as
k 11 k 12 k 14
k 22 k 23 0 k 26
k 33 k 34 0 k38
k 44 k 45 k 46 0
[K] = (112)
k 55 k 56 0
k 66 k 67 0
k 77 k 78
k 88
Friday 28 April 2000 – 14:43

In the skyline storage scheme, there is no equation for estimating

the number of elements below the skyline, and hence the storage
requirements. Nevertheless, the skyline storage requirements can
never be greater than the requirements of the symmetric banded
storage scheme. For the example at hand, the number of elements
under the skyline is 24. However, it should be noted that this
scheme requires an additional two auxiliary arrays. One for
storing the column heights of [K] ({a}) and one for storing the
location of the diagonal terms of [K] ({b}}. For the matrix of Eq.
112, these auxiliary arrays will contain the information shown in
Eq. 113. In computing the storage requirements of the skyline
storage scheme one needs to add the requirements of these two
arrays, i.e. 2n + 1 elements. Noting that these two arrays are of
type integer, their combined storage requirement in double
precision numbers can then be approximated as n. Therefore, for
the matrix of Eq. 112 the number of elements that need to be
stored is 24 + 8 = 32 .
From this example, one may conclude that for almost full
matrices the storage and addressing overhead associated with the
skyline storage scheme can outweigh its benefits. In general,
however, for matrices of the type encountered in finite-element
analyses, the skyline storage scheme and associated solvers can
be very effective. This is especially true if the node numbering of
the finite-element model is optimized so that the skyline size is
minimized. Node numbering optimization will be discussed in
the next section.
1 1
2 2
2 4
{a} = 4 n {b} = 6 n+1 (113)
2 10
5 12
2 17
6 19
25
Mesh The analyst using S-Frame is free to choose any convenient node
numbering scheme. However, once the finite-element model has
optimization
been defined, the S-Frame finite-element analysis engine
optimizes the internal node numbering so that the size of the
skyline is minimized. This node numbering optimization is
necessary for minimizing the memory requirements for storage
of the system matrices (stiffness, damping, mass).
The node-numbering optimization algorithm used in S-Frame is

due to Sloan (1986). This algorithm has been used on a broad
range of problems with various types of elements (e.g. problems
suggested by Everstine (1979)). The results suggest that it is
superior to both the reverse Cuthill-McKee (1969) and Gibbs-
King (1976) methods.
Internally, S-Frame works with three sets of node numbers:
1. User node numbers.
2. Internal node numbers.
3. Optimized node numbers.
To illustrate these node numbering schemes, let us consider the

simple finite-element model shown in Fig. 3.
‘
User 200 400 500 600 700 800 300
Internal 1 3 4 5 6 7 2
Optimized 1 2 3 4 5 6 7
Figure 3: Node numbering.

The user node numbers are arbitrary. They can start at any
number and they do not have to be sequentially numbered, that
is, different node numbering schemes may be used for different
parts of the structure. Of course a node number cannot be used
more than once, but gaps may exist in the numbers used. In
addition, all node numbers must be greater than zero. The
internal node numbers are assigned during the model definition.
In the example of Fig. 3, it is assumed that the analyst specified
Friday 28 April 2000 – 14:43

the nodes in the following order: 200, 300, 400, 500,..., 800.
Therefore, node 200 was specified first (internal node number 1),
node 300 was specified second (internal node number 2), node
400 was specified third (internal node number 3) and so on. The
optimized node numbers are computed by the node-numbering
optimization algorithm of S-Frame.
S-Frame uses an internal look-up table to convert from one

numbering scheme to another. All results and user messages use
the user node numbers, whereas all computations use the
optimized node numbers. Both the internal and the optimized
node numbers are transparent to the analyst.
Application of As noted earlier, once the system equations have been formed,
the system matrix (or matrices) has to be modified in order to
boundary
account for the boundary conditions. The system stiffness matrix
conditions will generally have a degeneracy equal in number to the rigid-
body modes of the structure. On application of adequate
boundary conditions, rigid-body modes and the degeneracy of
the stiffness matrix are removed and hence, the solution of the
system of equations becomes possible. The Boundary conditions
can be of two types: homogeneous and non-homogeneous. An
example of the first type is a constrained node, and of the second
type, a prescribed nodal displacement (foundation settlement).
Both types of boundary conditions can be applied using a
number of different schemes.
The scheme used in S-Frame is based on the concept of penalty

functions. A proof of this method is given in Zienkiewicz (1977).
Consider the system of equations
k 11 x 1 + k 12 x 2 + k 13 x 3 = r 1
k 21 x 1 + k 22 x 2 + k 23 x 3 = r 2 (114)
k 31 x 1 + k 32 x 2 + k 33 x 3 = r 3
and suppose that the second coefficient of the system nodal

displacement vector (x2) is prescribed the value of xp. Let us
further assume that the system of equations is symmetric and
that it is desirable to preserve this symmetry. We now select a
large number α which is six to twelve orders of magnitude larger
than the largest coefficient kij. We then add this coefficient to k22.
In addition, the right-hand-side of the second equation is

changed to α times the prescribed value, i.e. αxp. Following this
procedure, above system of equations becomes
k 11 x 1 + k 12 x 2 + k 13 x 3 = r 1
k21 x 1 + ( ( k 22 + α )x 2 ) + k23 x 3 = αx p (115)
k 31 x 1 + k 32 x 2 + k 33 x 3 = r 3
Since α is chosen to be much larger than all the equation

coefficients kij, for all practical purposes, the second equation is
equivalent to
αx 2 = αx p (116)
or
x2 = xp (117)
which is the desired result. This method of applying the

boundary conditions is very simple to implement, it conserves
symmetry, and allows the application of non-homogeneous
boundary conditions to be treated within load cases and/or load
combinations. The only disadvantage of this method is the
artificial introduction of very large diagonal terms which, in
certain cases, may introduce some ill-conditioning during the
solution phase, thereby resulting in possible loss of accuracy.
Therefore, if the finite-element model contains only homogenous
boundary condition, S-Frame uses a variation of the
aforementioned procedure which does not introduce the
unwanted ill-conditioning. Using this alternate procedure, after
applying the constraint
x2 = 0 (118)
Eqs. 114 become
k 11 x 1 + 0x 2 + k 13 x 3 = r 1
0x 1 + αx 2 + 0x 3 = 0 (119)
k 31 x 1 + 0x 2 + k 33 x 3 = r 3
Friday 28 April 2000 – 14:43

In S-Frame, the number α is set to the root-mean-square of the

diagonals of the system matrix, that is
k ii
α = -------------------
-. (120)
numdof
Constraints The only type of constraint supported by S-Frame is the coupling

of degrees of freedom. For this type of constraint, the analyst
(coupled degrees
must specify sets of nodes which have a certain degree of
of freedom) freedom coupled (i.e. these degrees of freedom are constrained
to have the same displacement, velocity and acceleration).
Consider, for example, the cantilever beam shown in Fig. 4.
Y
Ty direction
1 2 3 4 5 6 7 X
Rx direction
Figure 4: Constraints – coupled degrees of freedom.

In order to couple the translational y degrees of freedom of node
2, 3, 4 and the rotational x degrees of freedom of nodes 5, 6 the
user needs to specify the information given in Table 1.
In S-Frame, constraints can be specified among free degrees of

freedom only, i.e. degrees of freedom that have not been
restrained due to support conditions (boundary conditions
and/or support settlements). For instance, none of the six

degrees of freedom of node 1 in Fig. 4 (which is fully supported)
can be part of a constraint specification.
Node Coupled Degree

Numbers of Freedom
2 to 5 Ty
3 to 5 Ty
6 to 5 Rx
Table 1: Constraint – coupled degrees of freedom

The scheme used in S-Frame for application of constraints is as
follows. Consider the system of equations
k 11 x 1 + k 12 x 2 + k 13 x 3 = r 1
k 21 x 1 + k 22 x 2 + k 23 x 3 = r 2 (121)
k 31 x 1 + k 32 x 2 + k 33 x 3 = r 3
and suppose that the second degree of freedom (x2) is

constrained to be equal to the first degree of freedom (x1).
Furthermore, similar to the procedure for application of
boundary conditions, we assume that the system of equations is
symmetric and that it is desirable to preserve this symmetry. On
application of the constraint
x2 = x1 (122)
the system of equations 121 becomes
( k 11 + k 12 )x 1 + k 13 x 3 = r 1
( k 21 + k 22 )x 1 + k 23 x 3 = r 2 (123)
( k 31 + k 32 )x 1 + k 33 x 3 = r 3
Friday 28 April 2000 – 14:43

Since there are only two unknown (x1, x3), equations 123 are not
independent. If we take the second equation as the superfluous
one and add it to the first equation, then we are left with the
symmetric system of equations
( k 11 + k 12 + k 21 + k 22 )x 1 + ( k 13 + k 23 )x 3 = r 1 + r 2
(124)
( k 31 + k 32 )x 1 + k 33 x 3 = r 3
The main disadvantage of above procedure is that it does not

preserve the original size of the system of equations. Therefore,
in S-Frame the superfluous equations are not removed. Instead,
they are decoupled from the rest of the system of equations using
a technique similar to the one described for the application of
boundary conditions. For the example at hand, we first add the
2nd row to the first row so that equations 121 become
( k 11 + k 21 )x 1 + ( k 12 + k 22 )x 2 + ( k 13 + k 23 )x 3 = r 1 + r 2
k 21 x 1 + k 22 x 2 + k 23 x 3 = r 2 (125)
k 31 x 1 + k 32 x 2 + k 33 x 3 = r 3
and then we add the 2nd column to the first column to obtain
( k 11 + k 21 + k 12 + k 22 )x 1 + ( k 12 + k 22 )x 2 + ( k 13 + k 23 )x 3 = r 1 + r 2
( k 21 + k 22 )x 1 + k 22 x 2 + k 23 x 3 = r 2 (126)
( k 31 + k 32 )x 1 + k 32 x 2 + k 33 x 3 = r 3
All that remains is to decouple the second degree of freedom

from the rest of the equations. This is accomplished by zeroing
the 2nd row and column of the system of equations except for the
diagonal term, i.e.
( k 11 + k 21 + k 12 + k 22 )x 1 + 0x 2 + ( k 13 + k 23 )x 3 = r 1 + r 2
0x 1 + αx 2 + 0x 3 = α (127)
( k 31 + k 32 )x 1 + 0x 2 + k 33 x 3 = r 3
where α is again calculated according to Eq. 120. It can be easily

verified that above system of equations is equivalent to Eqs. 124.
Once the system of equations 127 is solved for x1, and x3, then
Eq. 122 is used to obtain the correct value for x2.
The analyst should be aware that the application of constraints

will normally increase the skyline storage requirements. Since
ultimately the constraints are applied using the optimized node
numbers, it is generally not possible for the analyst to estimate
this increase a priori. However, the following is a useful
guideline:
1. Coupling a lower degree of freedom (i) to a higher degree of

freedom (j) increases the memory requirements less than
coupling a higher degree of freedom (j) to a lower degree of
freedom (i).
2. The memory requirement increases as the numerical

difference between the coupled degrees of freedom
abs ( i – j )
increases.
Damping effects The most common form of viscous damping used in structural
dynamics is the so-called Rayleigh-type damping (Bathe (1982))
where the system damping matrix is expressed as a linear
combination of the system mass and stiffness matrices
[ C ]s = α [ M ] s + β [ K ] s (128)
The two damping parameters α and β may be determined by

specifying the damping ratio for two modes, say first and second,
to obtain
2ω 1 ω 2 ( ω 2 ξ 2 – ω 1 ξ1 )
α = ----------------------------------------------------
2 2
[ ω 2 – ω1 ]
(129)
2 ( ω 2 ξ 2 – ω1 ξ1 )
β = --------------------------------------
2 2
[ ω 2 – ω1 ]
The damping ratio in the other modes are then given by
α βω i
ξ i = -------- + --------- (130)
2ω i 2
Friday 28 April 2000 – 14:43

Mass-proportional and stiffness-proportional damping may also

be used separately. For mass-proportional damping β = 0.
Specifying ratio for the first mode only gives
α = 2ξ 1 ω 1 (131)
Therefore,
ξ1 ω1
ξi = ------------ (132)
ωi
which means that mass-proportional damping ratios decrease

with increase in mode frequency.
For stiffness-proportional damping α = 0. Again specifying the

damping ratio for the first mode only gives
2ξ 1
β = -------- (133)
ω1
Therefore,
ξ1 ωi
ξ i = ----------- (134)
ω1
In this case the damping ratios increase with increase of modal

frequency. In practice, it has been found that mass-proportional
damping can represent friction damping, whilst stiffness-
proportional damping can represent internal material damping.
Accepted values of modal damping ratios ξ i for typical forms of
construction vary from 0.01 for small diameter piping systems,
to 0.07 for bolted joint and reinforced concrete structures.
Damping effects may also be required for specific degrees of

freedom only. This is accomplished by attaching dashpots on
nodes in the desired directions. Dashpots are currently
supported for the Dynamic Analysis of S-Frame only.
The main advantage of Rayleigh damping is that the equations of

motion can be decoupled when using the normalized
eigenvectors (Eq. 33). In this case the generalized damping for
the ith mode is defined as
T
ci = { ψ } i [ C ] s { ψ }i (135)
Checking for non- Non-structural (unstable) degrees of freedom may arise in

S-Frame models as a result of nodal releases. The physical
structural degrees
interpretation of non-structural degrees of freedom is as follows:
of freedom
Non-structural degrees of freedom are degrees of freedom that
are not structurally connected to the rest of the structure.
They do not cause rigid body motion of the structure, however,

they do render the system stiffness matrix singular. Before the
solution phase of all types of analyses, S-Frame performs a check
to ensure that all degrees of freedom are connected to the
structure. If S-Frame finds a degree of freedom that is unstable
due to the imposition of nodal releases, it fixes the degree of
freedom by appropriate modifications to the stiffness matrix and
issues the CO0045 exception. If, however, S-Frame finds an
unstable degree of freedom which was not caused by the
imposition of nodal releases, it only issues the CO0046
exception.
Graph The S-Frame finite-element analysis engine often requires to

interpolate between user-specified data. For instance, in the case
interpolation
of a dynamic analysis, the engine must interpolate within
discrete data points of the user-specified force or acceleration vs.
time graph.
Friday 28 April 2000 – 14:43

S-Frame uses linear interpolation to compute the required value

of the dependent variable. Consider, for example, the data given
in Table 2 which is typical of the force vs. time graph required for
a dynamic analysis.
Time ( xi ) Force ( yi )
0.0 0.6
0.2 1.2
1.3 0.2
3.2 -0.3
5.6 -2.3
9.1 -1.1
Table 2: Discrete data.

During the time integration of the system equilibrium equations
46, the S-Frame finite-element analysis engine will require the
value of the force at times for which it has not been provided
explicitly by the user. For instance, if the value of the force ( y * ) at
time 1.0 (i.e. x * = 1.0 ) is required, S-Frame will interpolate using
the following equation
y2 – y1
y * = y 1 + ---------------- ( x * – x 1 ) (136)
x2 – x1
or
0.2 – 1.2
y * = 1.2 + -------------------- ( 1.0 – 0.2 ) = 0.47272 (137)
1.3 – 0.2
Lumped mass VS. The consistent mass matrix is based on the approximation for the
kinetic energy. That is, while the approximation for the strain
consistent mass
energy leads to the stiffness matrix, the use of the same
formulations polynomial functions (or shape functions) in the approximation
of the kinetic energy, leads to the consistent mass matrix. As a
result one finds some negative off diagonal elements in the
consistent mass matrix. The consistent mass matrix is a more
accurate representation of the inertial properties and it further
leads to natural frequencies that are bounded below by the exact
natural frequencies. That is, the computed natural frequencies

always lie above the exact values and by using sufficient number
of elements one can obtain very accurate results.
Lumped mass representation, which predates the introduction of

the consistent mass matrix, is based on equivalent mass
representation in the gravitational field. Thus, in the case of the
beam element one may choose to place half the mass of the beam
at each node of the element. However, such representation
overlooks the inertial effects due to flexure and in rigid-body
rotations of the beam it will over estimate the rotational inertia.
Accordingly, the lumped mass representation is not as accurate
and can lead to natural frequencies that may be higher or lower
than the exact ones. However, for large problems lumped mass
representation leads to considerable savings in computing since
it gives rise to a diagonal mass matrix. The current release of
S-Frame supports the consistent mass formulation only. In a
future release the analyst will be given the option to use either
consistent or lumped mass formulation.
Number of The cubic polynomial used to develop the element matrices for
the beam element is an exact solution of the homogeneous part
elements required
of the beam static equilibrium equations. Accordingly, the
stiffness matrix is exact and often one need only use one beam
element per span in static problems. In non-static problems
(vibration, dynamic, buckling) the equilibrium equations are
different and the cubic polynomial is no longer an exact solution.
For instance, the solution of the free-vibration problem is
frequency dependent. Accordingly, for non-static analyses one
must use several elements to obtain accurate results. This is
especially true in the following cases:
1. Vibration: higher frequencies are of interest.
2. Buckling: higher buckling modes are of interest.
3. Dynamic: forcing function includes high frequencies.
4. Response Spectrum: spectrum includes high frequencies.
In cases where the model includes other type of elements (e.g.

shell elements) above discussion is also true for static analysis.
Friday 28 April 2000 – 14:43

S-Frame Theory Manual Chapter 5 : Element library : 62
5 Element library
S-Frame’s element library contains the following elements:
1. A beam element (with axial, bending, and torsional

stiffness/inertia).
2. A thin flat shell element (with membrane and bending

stiffness/inertia).
3. A truss element (with axial stiffness and translational

inertia).
4. A two-noded spring (with axial and torsional stiffness).
5. Tension-only members
6. Compression-only members
7. Spring element.
8. Hook element.
9. Gap element.
63 : Chapter 5 : Element library S-Frame Theory Manual
Above elements are used in linear and nonlinear analysis. In the

case where a an element type is not appropriate for a certain
analysis type (e.g. a cable element in a linear static analysis), the
solver performs an automatic element and spring conversions
according to Table 3 and Table 4.
Stressed Vibration
Nonlinear Static
Free Vibration
Dynamic
Buckling
P-Delta
R.S.A.
Static
Element
Type
Description
0 Non-active 0 0 0 0 0 0 0 0
1 Beam 1 1 1 1 1 1 1 1
2 Truss 2 2 2 2 2 2 2 2
3 Cable 2 2 2 2 2 2 2 3
4 2-noded spring 4 4 4 4 4 4 4 4
5 2-noded torsion 5 5 5 5 5 5 5 5
spring
6 2-noded nonlinear 4 4 4 4 4 4 4 6
spring
7 2-noded nonlinear 5 5 5 5 5 5 5 7
torsion spring
8 Tension-only truss 2 2 2 2 2 2 2 8
9 Compression-only 2 2 2 2 2 2 2 9
truss
Table 3: Automatic solver element conversions
Friday 28 April 2000 – 14:43

Stressed Vibration
Nonlinear Static
Free Vibration
Dynamic
Buckling
P-Delta
R.S.A.
Static
Spring
Type
Description
0 Compression only 2 2 2 2 2 2 2 0
if motion -ve
1 Compression only 2 2 2 2 2 2 2 1
if motion +ve
2 Linear spring 2 2 2 2 2 2 2 2
3 Nonlinear spring 2 2 2 2 2 2 2 3
Table 4: Automatic solver spring conversions

The following sections describe these elements in more detail.
3D Beam The 3D-beam element is a uniaxial element with tension,

compression, torsion, and bending capabilities. It has six degrees
of freedom at each node: three translations in the element x, y, z
directions and rotations about the element x, y, z axes.
Stiffness, geometric The finite element formulation of the 3D-beam element uses a
stiffness, and mass third order polynomial for bending and linear polynomials for
matrices tension and torsion (Przemieniecki (1968), Cook (1989), Weaver
(1980)). The stiffness, geometric stiffness and consistent mass
matrices are given by Eqs. 138, 147 and 148 respectively.
AE
------
L
1
0 kz
1
0 0 ky
GJ
0 0 0 ------
L
4 5
0 0 ky 0 ky
3 5
0 kz 0 0 0 kz
[ K ]e = (138)
A E-
–-------- AE
0 0 0 0 0 ------
L L
2 4 1
0 kz 0 0 0 kz 0 kz
2 3 1
0 0 ky 0 k y 0 0 0 ky
– GJ GJ
0 0 0 --------- 0 0 0 0 0 ------
L L
4 6 3 5
0 0 ky 0 k y 0 0 0 kz 0 k y
3 6 4 5
0 kz 0 0 0 kz 0 ky 0 0 0 k z
Friday 28 April 2000 – 14:43

where
s
Ai shear area normal to direction i.
L element length.
G shear modulus.
J torsional moment of inertia.
1 12EIi
k i = -----------------------
3
- (139)
L ( 1 + Γi )
2 – 12 EIi
k i = -----------------------
3
- (140)
L ( 1 + Γi )
3 6EIi
k i = -----------------------
2
- (141)
L ( 1 + Γi )
4 – 6 EI i
k i = -----------------------
2
- (142)
L ( 1 + Γi )
5 ( 4 + Γ i )EI i
k i = ------------------------- (143)
L ( 1 + Γi )
6 ( 2 – Γ i )EI i
k i = ------------------------ (144)
L ( 1 + Γi )
12EIz
Γ y = --------------
2 s
(145)
GL A z
12EI y
Γ z = --------------
2 s
- (146)
GL A y
0
0 36
0 0 36
0 0 0 0
2
0 0 – 3L 0 4L
Pe 2
[ K g ] e = ----- 0 3L 0 0 0 4L (147)
30 0 0 0 0 0 0 0
0 – 36 0 0 0 – 3L 0 36
0 0 – 36 0 3L 0 0 0 36
0 0 0 0 0 0 0 0 0 0
2 2
0 0 – 3L 0 – L 0 0 0 3L 0 4L
2 2
0 3L 0 0 0 – L 0 – 3L 0 0 0 4L
where
P axial load (positive for a tensile load).
140
0 156
0 0 156
140J
-----------
0 0 0
A
2
0 0 – 22L 0 4L
me 2
[ M ] e = -------- 0 22L 0 0 0 4L (148)
420 70 0 0 0 0 0 140
0 54 0 0 0 13L 0 156
0 0 54 0 – 13L 0 0 0 156
70J 140J
0 0 0 -------- 0 0 0 0 0 -----------
A A
2 2
0 0 13L 0 –3L 0 0 0 22L 0 4L
2 2
0 –13L 0 0 0 – 3L 0 – 22L 0 0 0 4L
where
me element mass.
For the nonlinear beam element, the 2nd order terms are also
included in the tangent stiffness matrix.
Friday 28 April 2000 – 14:43

3D Thin Shell The triangular and quadrilateral thin shell elements have both
out-of-plane (bending) and in-plane (membrane) capabilities.
Both in-plane and out-of-plane loads are permitted. These
elements have six degrees of freedom at each node: three
translations in the element x, y, z directions and rotations about
the element x, y, z axes. For nonlinear static analysis, the 2nd
order terms are also included in the tangent stiffness matrix of
the thin shell family of elements.
Stiffness and mass In the element coordinate system (with x, y in the plane of the
matrices element and z normal to it) the in-plane and out-of-plane actions
are decoupled and one may derive the stiffness and mass
matrices, for in-plane and out-of-plane actions, independently.
The finite-element formulation of these thin shell elements uses a
linear polynomial for the in-plane action, and a third order
polynomial for the out-of-plane actions. Unlike in the case of the
beam element, the element matrices (stiffness and consistent
mass) for the shell elements cannot be evaluated explicitly and
they require the use of numerical integration algorithms. For
more details on these shell elements, the user is referred to
Bazeley et. al. (1965). A Softek internal report on the thin shell
element formulation is available.
Plate elements The triangular and quadrilateral thin plate elements have out-of-
plane (bending) stiffness only. Only out-of-plane loads are
permitted. These elements have three degrees of freedom at each
node. With x, y in the plane of the element and z normal to it, the
nodal degrees of freedom are: translation in the z direction, and
rotations about the x and y directions.
Membrane elements The triangular and quadrilateral thin membrane elements have
in-plane (membrane) stiffness only. Only in-plane loads are
permitted. These elements have two degrees of freedom at each
node. With x, y in the plane of the element and z normal to it, the
nodal degrees of freedom are: translations in the x and y
directions.
Aspect ratio It is important that the aspect ratio of these thin shell, plate and
membrane elements (defined as the ratio of the longest side to
the shortest side of an element) be lower than 4. Highly distorted
elements are known to yield poor results.
3D Truss The 3D-truss element is a uniaxial element with tension

capabilities only. It has one degree of freedom at each node: a
translation in the element x direction. For the nonlinear truss
and cable elements, the 2nd order terms are also included in the
tangent stiffness matrix. A Softek internal report on the
nonlinear truss and cable element formulation is available.
Stiffness and mass The finite element formulation of the 3D-truss element uses a
matrices linear polynomial for tension (Cook (1989), Weaver (1980)). The
stiffness and consistent mass matrices are given by Eqs. 149 and
150 respectively.
AE
------
L
0 0
0 0 0
0 0 0 0
0 0 0 0 0
[ K ]e = 0 0 0 0 0 0 (149)
–A E AE
--------- 0 0 0 0 0 ------
L L
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
Friday 28 April 2000 – 14:43

where
L element length.
140
0 140
0 0 140
0 0 0 0
0 0 0 0 0
me 0 0 0 0 0 0
[ M ] e = -------- (150)
420 70 0 0 0 0 0 140
0 70 0 0 0 0 0 140
0 0 70 0 0 0 0 0 140
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
where
me element mass.
Two-noded spring The two-noded spring element has tension and torsional
stiffness. It has two degrees of freedom at each node: one
translation in the element x direction and a rotation about the
element x axis.
Stiffness matrix The stiffness matrix is given by
ka
0 0
0 0 0
0 0 0 kr
0 0 0 0 0
[ K ]e = 0 0 0 0 0 0 (151)
–ka 0 0 0 0 0 k a
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 – kr 0 0 0 0 0 k r
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
where
ka axial stiffness of linear spring.
kr torsional stiffness of linear spring.
For the nonlinear 2-noded spring, k a and k r are not constant.

They are functions of the spring deformation (axial and
rotational respectively) provided in the form of a table.
Tension-only Tension-only members are truss elements with a stiffness matrix

give by Eq. 151
members
where
if member is in compression k a = 0
otherwise k a = k a
kr = 0
Therefore, as their name implies, tension-only members have

stiffness only when they are in tension.
Compression-only Compression-only members are truss elements with a stiffness

matrix give by Eq. 151
members
Friday 28 April 2000 – 14:43

where
if member is in tension k a = 0
otherwise ka = k a
kr = 0
Therefore, as their name implies, compression-only members

have stiffness only when they are in compression.
Spring element The spring element can have tension and torsional stiffness. It
has two degrees of freedom: one translation in the element x
direction and a rotation about the element x axis.
Stiffness matrix The stiffness matrix is given by
ka
0 0
[ K ]e = 0 0 0 (152)
0 0 0 kr
0 0 0 0 0
0 0 0 0 0 0
where
ka axial stiffness of linear spring.
kr torsional stiffness of linear spring.
For nonlinear static analysis, k a and k r are not constant. They are
functions of the spring deformation (axial and rotational
respectively) provided in the form of a table.
Hook element The hook element is a spring element which engages after it
undergoes a compression greater or equal to the user specified
GAP. Its stiffness is given by Eq. 152
where
if spring compression < GAP ka = 0
otherwise k a = k a
kr = 0
Gap element The gap element is a spring element which engages after it
undergoes an elongation greater or equal to the user specified
GAP. Its stiffness is given by Eq. 152
where
if spring elongation < GAP k a = 0
otherwise k a = k a
kr = 0
Friday 28 April 2000 – 14:43

S-Frame Theory Manual Chapter 6 : Unit conversion factors : 74
6 Unit conversion factors

Unit conversions in S-Frame are performed during the model
definition and throughout the analysis. The following six tables
give the unit conversion factors used by S-Frame.
Second
Second (s) 1.0
Minute (min) 60.0
Hour (h) 3600.0
Table 5: Unit conversion factors for time
Pound
Pound (lbf) 1.0
Ounce (US oz) 0.0625
Kip (kip) 1000.0
Ton (tnf) 2000.0
Kilogram (kgf) 2.20463414
Tonne (tf) 2204.63414
Dyne (dynf) 2.248089237E-6
Newton (N) 0.2248089237
Kilonewton (kN) 224.8089237
Meganewton (mN) 224808.9237
Table 6: Unit conversion factors for force
Inch
Inch (in) 1.0
Foot (ft) 12.0
Yard (yd) 36.0
Table 7: Unit conversion factors for length
75 : Chapter 6 : Unit conversion factors S-Frame Theory Manual
Inch
Millimeter (mm) 0.03937007874
Centimeter (cm) 0.3937007874
Meter (m) 39.37007874
Table 7: Unit conversion factors for length
Mug
Mug (mg) 1.0
Pound (lbm) 0.00259041
Kip (kipm) 2.59041
Ton (tnm) 5.18082
Tonne (tm) 5.71081
Slug (-) 0.0833333
Gram (g) 5.71081E-6
Kilogram (kg) 5.71081E-3
Table 8: Unit conversion factors for mass
Radian
Radian (rad) 1.0
Degree (o) 0.0174532925
Cycle (cycle) 6.283185307
Table 9: Unit conversion factors for angle
oF
oF 1.0
o
C 1.8
Table 10: Unit conversion factors for temperature
Friday 28 April 2000 – 14:43

S-Frame Theory Manual Chapter 6 : Unit conversion factors : 76
It should be noted that the S-Frame finite-element analysis

engine converts all quantities to the International System of
Units, abbreviated SI (from the French, Système International
D’Unités), prior to the solution phase. Table 11 gives the basic SI
units used by S-Frame.
Quantity Basic SI Unit

Length meter (m)
Time second (s)
Mass kilogram (kg)
Force newton (N)
Table 11: SI system of units
77 : Chapter 6 : Unit conversion factors S-Frame Theory Manual
Friday 28 April 2000 – 14:43

S-Frame Theory Manual Chapter 7 : Memory requirements : 78
7 Memory requirements
The following sections describe the memory requirements (both
Random Access Memory (RAM) and hard disk storage) of the
S-Frame solution engine.
Dynamic memory The following expressions give an approximation to the number

of bytes of RAM the S-Frame solution engine must dynamically
allocation
allocate during execution.
requirements
(RAM)
Linear static: numdof ( 4 numdof + 32 ) (153)
Nonlinear static: numdof ( 4 numdof + 40 ) (154)

2
Free vibration: numdof ( 8 numdof + 8 numitrvec + 24 ) + 24 numitrvec (155)
2
Stressed free-vibration: numdof ( 8 numdof + 8 numitrvec + 32 ) + 24 numitrvec (156)
Linear dynamic: numdof ( 16 numdof + 120 ) (157)
P-delta: numdof ( 8 numdof + 32 ) (158)

2
Linear buckling: numdof ( 12 numdof + 8 numitrvec + 40 ) + 24 numitrvec (159)
2
numdof ( 8 numdof + 8 numitrvec + 24 ) + 24 numitrvec +
Response spectrum: (160)
2 ( 96 num1Delm + 296 num2Delm )
Storage The following expressions give an approximation to the number

of bytes of permanent storage the S-Frame solution engine needs
requirements
during execution.
Input files: 2 ( 62 numnod + 62 num1Delm + 31 num2Delm ) (161)

N ( 68 numnod + 72 num1Delm + 164 num2Delm + 8073 ) +
Output files: 64 ( numnod numnatfrq ) + (162)
296 ( numRSAldcas num2Delm numnatfrq )
where
N = Max ( numldcas + numldcmb, numtimstp ) (163)
79 : Chapter 7 : Memory requirements S-Frame Theory Manual
Temporary files: Linear and nonlinear static:

( 7512 + 196 numldcas ) ( num1Delm + num2Delm ) + 48 numnod +
(164)
numdof ( 4 numdof + 8 numldcas )
Free vibration:
2
7512 ( num1Delm + num2Delm ) + 48 numnod +8 numdof (165)
Stressed free-vibration:
(166)
Linear dynamic:
(167)
P-delta:
(168)
Linear buckling:
(169)
Response spectrum:
2
7512 ( num1Delm + num2Delm ) + 48 numnod +8 numdof (170)
In order to minimize temporary storage when performing RSA

analysis, it is highly recommended that the user makes the RSA
load cases the first ones.
It should be noted that all temporary files are removed by the

S-Frame solution engine after the completion of the analysis.
Friday 28 April 2000 – 14:43

S-Frame Theory Manual Chapter 8 : References : 80
8 References
Anton, H., Elementary Linear Algebra, John Wiley & Sons,
Toronto, 1977.
Bathe, K-J, Finite Element Procedures in Engineering Analysis,

Prentice-Hall Inc., Englewood Cliffs, New Jersey, 1982.
Bazeley, G.P., Cheung, Y.K., Irons, B.M., Zienkiewicz, O.C.,

Triangular Elements in Plane Bending – Conforming and Non-
conforming Solutions, Proceedings 1st on Matrix Methods in
Structural Mechanics, Ohio, pp. 547-576, 1965.
Cathill, E., McKee, J., Reducing the Bandwidth of Sparse

Symmetric Matrices, Proceedings ACM National Conference,
Association of Computing Machinery, New York, 1969.
Clough, R.W., Penzien, J., Dynamics of Structures, McGraw-Hill,

Toronto, 1975.
Cook, R.D., Malkus, D.S., Plesha, M.E., Concepts and

Applications of Finite Element Analysis, John Wiley & Sons,
Toronto, 1989.
Craig, R.R., Structural Dynamics, An Introduction to Computer

Methods, John Wiley & Sons, Toronto, 1981.
Everstine, G.C., A Comparison of three Resequencing Algorithms

for the Reduction of Matrix Profile and Wavefront, International
Journal of Numerical Methods in Engineering, 14, 837-853, 1979.
Gerald, C.F., Applied Numerical Analysis, 2nd ed., Addison-

Wesley, Don Mills, Ontario, 1980.
Gibbs, N.E., A Hybrid Profile Reduction Algorithm, ACM Trans.

Math. Software, 2, 378-387, 1976.
Jennings, A., Matrix Computation for Engineers and Scientists,

John Wiley & Sons, 1977.
Petyt, M., Introduction to Finite Element Vibration Analysis,

Cambridge University Press, New York, 1990.
81 : Chapter 8 : References S-Frame Theory Manual
Press, W.H., Flannery, B.P., Teukolsky, S.A., Veuerling, W.T.,

Numerical Recipes in FORTRAN, Cambridge University Press,
Cambridge, New York, 1990.
Przemieniecki, J.S., Theory of Matrix Structural Analysis,

McGraw-Hill, Toronto, 1968.
Sloan, S.W., An Algorithm for Profile and Wavefront Reduction of

Sparse Matrices, International Journal of Numerical Methods in
Engineering, 23, 239-251, 1986.
Stasa, F.L., Applied Finite Element Analysis for Engineers, Holt,

Rinehart and Winston, Toronto, 1985.
Tabarrok, B., Stylianou, M.C., Comparison of Integration Schemes

for Vibration of Structures, CSME Proceedings of the 10th
Canadian Congress of Applied Mechanics, D67-68, 1985.
Weaver, W., Gere, J.M., Matrix Analysis of Framed Structures, D.

Van Nostrand Company, Toronto, 1980.
Weaver, W., Johnson, P., Structural Dynamics by Finite Elements,

Prentice-Hall Inc., Englewood Cliffs, New Jersey, 1987.
Zienkiewicz, O.C., The Finite Element Method, McGraw-Hill

(UK), London, 1977, pp. 204-207.
Friday 28 April 2000 – 14:43

S-Frame Theory Manual

Uploaded by

Copyright:

Available Formats

S-Frame Theory Manual

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

S-Frame Theory Manual

Uploaded by

Copyright:

Available Formats

What analysis procedures are described in the document?

What analysis procedures are described in the document?

What elements are included in the element library?

What elements are included in the element library?

S-Frame

Friday 28 April 2000 – 14:40

Unit conversion factors . . . . . . . . . . . . 74

Friday 28 April 2000 – 14:40

Friday 28 April 2000 – 14:40

Friday 28 April 2000 – 14:43

Ri RSA Base Shear in ithdirection.

Si Scaling factor in ithdirection.

[A] a square matrix.

[ K ]s a shifted system stiffness matrix.

Friday 28 April 2000 – 14:43

{ η }s generalized system displacement vector.

numdir number of ground motion directions.

Friday 28 April 2000 – 14:43

then their contribution to the behavior of the overall structure

to yield the system equations

where the internal load vector, for linear analysis, is defined as

linear static analysis

linear dynamic analysis

nonlinear static analysis

Linear static The purpose of the static analysis is to determine the

1. Stiffness effects and applied loads do not depend on time.

2. Inertial and damping effects are ignored.

3. Static acceleration fields, such as gravity, may be included

4. Time independent loads, displacements, pressures and

where the system load vector { R } s includes the contributions of:

1. Applied nodal loads.

2. Loading due to static acceleration fields (such as gravity).

3. Element thermal/pressure loads.

The system stiffness matrix will have a degeneracy equal in

and a unique solution { X } s may be obtained. Once the nodal

Friday 28 April 2000 – 14:43

Implementation S-Frame uses a skyline implementation (see Skyline storage

The stiffness matrix factorization, a computationally intensive

During matrix factorization, the skyline solver performs the

1. If a pivot is less than zero (system matrix is negative definite),

2. If a pivot is identically zero (system matrix is positive semi-

3. If a pivot is greater than zero but less than tstval, then it

Free-vibration The purpose of the free-vibration analysis is to determine the

1. Stiffness and inertial effects of the structure are time

2. Damping effects are ignored.

3. No forces, displacements, pressures or temperature effects

In a free-vibration analysis the motion is simple harmonic and

The acceleration vector may be obtained by taking the second

Substituting Eqs. 13 and 14 into Eq. 12 yields the following

Friday 28 April 2000 – 14:43

which may be re-written as

Hence, a non-trivial solution is possible only when the

which is the basic statement for free vibration problems. The

for the numdof eigenpairs

Subspace iteration is a very effective method of determining the

1. The number of eigenpairs (numnatfrq) to be extracted is set

numdof = 6 numnod – numbc (21)

2. Select a starting matrix [ A ]1 of dimension (numdof x

numitrvec = Min ( 2 numnatfrq, numnatfrq + 8 ) (22)