NONLINEAR ANALYSIS OF DYNAMIC STABILITY
OF ELASTIC SHELLS OF REVOLUTION
By
MARCUS GEORGE HENDRICKS
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA IN PARTIAL
FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1974
TO MY WIFE,
GENEVIEVE,
AND MY DAUGHTER,
CARMEN
ACKNOWLEDGMENTS
The author wishes to express his sincere appreci-
ation to Dr. S. Y. Lu, chairman of his supervisory com-
mittee, for his continuous guidance and encouragement
throughout the entire period of this study. He also
wishes to express his gratitude to Dr. I. K. Ebcioglu,
Dr. U. H. Kurzweg, Dr. R. L. Sierakowski, and Dr. M. W.
Self for their helpful discussions with the author and
many valuable suggestions.
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS ................................. iii
LIST OF FIGURES................................. vi
NOMENCLATURE .................................... vii
ABSTRACT ........................................ x
CHAPTER
1 INTRODUCTION ............................ 1
2 THEORETICAL DEVELOPMENT OF
NONLINEAR SHELL EQUATIONS ............... 14
2.1 Shell Coordinates.................. 14
2.2 Kinematic Relations................ 15
2.3 Constitutive Relations............. 17
2.4 Kirchhoff's Hypothesis............. 18
2.5 Variational Principle.............. 20
2.5.1 Hamilton's Principle ........ 20
2.5.2 Equations of Motion......... 21
2.5.3 Boundary conditions ......... 26
2.6 Synthesis of Equations............. 26
2.6.1 Shells of Revolution ........ 27
2.6.2 Conical Shells .............. 27
3 DYNAMIC STABILITY EQUATIONS
OF CONICAL SHELLS ....................... 35
4 METHOD OF SOLUTION...................... 39
4.1 Circumferential Modal Analysis ..... 39
4.2 Elimination of
Meridional Derivatives............. 42
4.2.1 Subdomain Method............ 45
4.2.2 Method of Satisfying
Boundary Conditions ......... 47
TABLE OF CONTENTS (Continued)
Page
5 NUMERICAL ANALYSIS FOR CONICAL SHELLS..... 52
5.1 Functional Description of the
Computer Program..................... 52
5.2 Checking the Computer Program........ 54
5.3 Dynamic Stability Results............ 55
5.3.1 Dynamic Modal Preference...... 56
5.3.2 Linearized Dynamic Stability.. 58
5.3.3 Nonlinear Dynamic Stability... 60
6 RESULTS AND CONCLUSIONS................... 69
APPENDIX I SHELL SPECIMEN AND BOUNDARY
CONDITIONS ............................ 71
APPENDIX II CONVERGENCE .......................... 74
APPENDIX III FLOW DIAGRAM OF THE
COMPUTER PROGRAM .................... 77
APPENDIX IV COMPUTER PROGRAM SOURCE
LISTING .............................. 81
BIBLIOGRAPHY ....................................... 126
BIOGRAPHICAL SKETCH ................................ 131
LIST OF FIGURES
Figure Page
1 Development Paths in Constructing
Solution Methods for Problems of
Dynamic Stability of Shells............. 4
2 Shell Element........................... 25
3 Conical Shell Coordinates............... 29
4 Fourier Wave Number Dependence.......... 57
5 Flexural Mode Preference During
Dynamic Instability..................... 59
6 Linearized Dynamic Stability for
High Axial Load ......................... 61
7 Linearized Dynamic Stability for
Intermediate Axial Load ................. 62
8 Nonlinear Dynamic Stability for
Heaviside Pressure Loading.............. 65
9 Nonlinear Dynamic Stability for
Intermediate Axial Load................. 66
10 Dynamic Stability Interaction Curves.... 68
11 Convergence Properties of Solution ...... 76
12 Flow Diagram of the Computer Program.... 80
NOMENCLATURE
A,B
[B]
D
E
G
K
L
Mr ,MB,MaB
N
N ,N ,N
P
P
Q ,Q
R
Ro
T
U,V,W
W
Coefficients of the first fundamental form
Boundary condition matrix
Flexural stiffness
Young's modulus
Shear modulus
Membrane stiffness
Lagrangian Function
Stress couple components
Stress couple components at edges
Number of terms in power series
Stress resultants
Stress resultants at edges
Axial load
Potential energy
Shear stress resultants
Position Vector
Radius of latitude circle
Principal radii of curvature
Kinetic energy
Displacements of arbitrary point in the shell
Displacement function for the buckling mode
vii
V
an ,b ,c n
a a ,an
e ,e ,e
f,g,h
h
n
j
m
n
a 6n
p ,p ,p
Pn
q
r
s
t
u,v,w
x3
z
a,B
Y
C
Region
Temporal coefficients
Coefficients of thermal expansion
Strain components
Differential operators
Shell thickness
Normal unit vector
Number of essential boundary constraints
Number of time dependent unknowns
Circumferential wave number
Surface loads
Fourier component of normal pressure
Transformed coordinate
Radius vector
Meridional conical coordinate
Nondimensionalized meridional conical
coordinate
Time
Middle surface displacements
Nondimensionalized middle surface
displacements
Coordinate of axis of revolution
Distance along i
Curvalinear coordinates or conical angles
Strain tensor
Error
viii
Circumferential conical coordinate
Eigenvalue
Poisson's ratio
Density
Stress tensor
Nondimensional time
Spatial coefficient
Transformed circumferential conical
coordinate
Additional
Notations
[ ]
[ ]
( ),
(') ,( )
L J
{ }
Number in brackets refers to reference
numbers in the Bibliography
Letter in brackets denotes matrix
Comma denotes partial differentiation
Denotes partial differentiations with
respect to time
Row vector
Column vector
Abstract of Dissertation Presented to the
Graduate Council of the University of Florida in Partial Fulfillment
of the Requirements for the Degree of Doctor of Philosophy
NONLINEAR ANALYSIS OF DYNAMIC STABILITY
OF ELASTIC SHELLS OF REVOLUTION
By
Marcus George Hendricks
December, 1974
Chairman: Dr. S. Y. Lu
Major Department: Engineering Sciences
It is the purpose of this dissertation to investi-
gate the dynamic stability of elastic shells of revolu-
tion. Two specific areas of this broad field are treated
in detail. First, this analytical study generates pre-
viously unavailable interaction relations for combined
dynamic loadings as they interact to cause passage from a
dynamically stable state to other states. Secondly, the
concept of linearized dynamic stability is extended to
include geometrically nonlinear effects. One role of these
effects is to allow the possibility of autoparametric ex-
citation of preferred flexural modes by the driven modes.
The question of whether such excitation can occur during
the primary dynamic instability motion sufficiently to af-
fect the magnitude of critical dynamic loads is studied.
A modified version of the subdomain method in combination
with circumferential modal analysis is developed for the
solution method. A computer program is constructed to
obtain numerical results and the dynamic stability charac-
teristics of a conical frustum are studied under a variety
of combined dynamic loadings. Interaction curves for
static stability, linearized dynamic stability, and non-
linear dynamic stability are generated.
For the particular loading conditions studied, the
results indicate that dynamic instability will occur at
loads below critical static loads. This reduction in crit-
ical dynamic loads was shown to be the result of both the
dynamic load factor effect and of autoparametric excitation.
The interaction curves which are generated illustrate these
effects quantitatively under conditions of combined dynamic
loading. A criterion for dynamic buckling is established
based on meridional mode shape changes. This ability to
detect sudden jumps in the meridional profile helps to verify
instability detected by divergence on displacement time
history curves and provides additional information about
the poststable state.
CHAPTER 1
INTRODUCTION
A dynamic stability problem can be defined as a
problem analyzed by Newton's equations of motion or any
equivalent method [1]. Dynamic stability problems for
continuous systems are nearly always governed by nonlinear
partial differential equations [2]. The solution to a
problem of dynamic stability in thin-walled structures may
be stated as the determination of what load-time combina-
tions cause the displacements of points in the material
body to increase sufficiently to either interfere with
operational specifications or cause a breakdown of the
structure. Methods for obtaining such solutions may be
categorized into three fundamental areas. The most obvi-
ous and usually most difficult method is the direct inte-
gration of the governing equations for sufficiently long
periods to allow direct observation of the ensuing motions.
A second technique is the investigation of the stability
of equilibrium points. The question here is whether a
slight perturbation of a dynamical system from an equilib-
rium state will produce motion confined to the neighborhood
of the equilibrium point or a motion tending to leave
that neighborhood [3]. This method is usually based on
either perturbation or characteristic equation techniques.
A third category of dynamic stability methods is on a
higher level of abstraction than the previous two and
relies on the existence and properties of variously de-
fined functionals. Energy functionals [4] and Liapunov
functionals [5] are two examples.
Definitions of dynamic stability may be categorized
by methods of approach to such problems. In situations
where the direct integration approach is employed, dynamic
instability may be defined to exist when some character-
istic parameter, such as deflection, stress, or strain,
becomes unbounded in one or more parts of the body. When
the method of small oscillations about equilibrium points
is used, instability may be defined as unbounded growth
of perturbations with time or by the existence of one or
more positive real parts of roots of a characteristic
equation. In problems where the solution method is of
the third above mentioned category a sufficient condition
for stability is defined according to the definite form
of the functional employed. In the case of energy func-
tionals, a boundary between stable and unstable regions
is defined by the vanishing of the second variation of
the energy functional. When Liapunov functionals are
employed, stability is defined when the functional is
positive definite and when its derivative is negative
definite.
The particular problem of interest in this study
is the dynamic buckling of shell structures, which rep-
resents a specific application of the more general theory
of dynamic stability principles. Just as the definitions
of dynamic stability are related to the solution methods,
the definitions of dynamic buckling of shells are related
to solution methods as well as the specifics of the phys-
ical problem such as geometry, material properties, or
loading rates.
Solution methods, falling within the first of the
three aforementioned categories, are direct temporal in-
tegration of equations of motion. They may be classified
into five procedural divisions. These are finite differ-
ence, finite element, Hamilton's principle, energy methods,
and methods of weighted residuals. Figure 1 illustrates
the procedural flow and interrelations between these meth-
ods for the study of dynamic stability of elastic shells
and is based on the most prominent articles dealing with
this subject.
First, let us consider recent significant contri-
butions falling within the finite difference method divi-
sion. Witmer et al. [6, 7, 8] have developed during the
1_
Theory of Theory of Newton's Laws
Elasticity Surfaces or
Hamilton's Variati
Principle
Coupled Set of
Nonlinear Partial
Differential
Equations of
Equilibrium
---- tL
Modal Analysis
to Eliminate
Circumferential
Coordinate
I
- I
Coupled Set of
Nonlinear Partial
Differential
Expressions of
Internal Energy
--- *-----
Finite Difference
in Meridional and
Circumferential
Coordinates
Assumed Solution Finite
in Meridional in Meri
Coordinate Coordir
-it
Difference
dional
tate
Kinematic
Relations
onal of a
Finite
Element
Assemblage of
Elements in
Dynamic
Equilibrium
Coupled Nonlinear System(s)
of Ordinary Differential
Equations or of
Integro-Differential
Equations
Numerical Integration i
Shallow shell theory is often introduced here
Bifurcation, if enforced, is usually introduced here
Figure 1 .
Development Paths in Constructing Solution
Methods for Problems of Dynamic Stability
of Shells
L
past decade a shell dynamics computer program known as
PETROS which utilizes spatial finite difference methods
in both the circumferential and meridional directions.
Utilizing PETROS the authors have studied a number of
explosive forming problems and have shown remarkable
agreement with experimental results. However, consider-
ing that after a ten-year development period, computer
run times are still measured in hours for moderately com-
plex problems, the practical usefulness of PETROS in study-
ing shell dynamic stability is somewhat limited. Using
shallow shell theory, Longhitano [9] developed numerical
solutions for the nonlinear dynamics of hemispherical shells
also utilizing finite difference solutions in two directions.
He showed that suddenly applied symmetic pressures may
couple the breathing mode into more preferred axisymmetric
or asymmetric vibratory modes causing instability at approx-
imately one half static buckling pressures. Again, however,
the two-dimensional discrete type of spatial treatment re-
sults in unusually long computer run times thus limiting
the number of cases treated.
Another class of finite difference solutions to
shell dynamic stability problems have appeared in the lit-
erature which require only approximately three fourths of
the computer run times as the previously cited two-dimen-
sional discrete methods and appear to give comparable
results. Such solutions based on elimination of the
circumferential coordinate via modal analysis prior to
finite differencing in the meridional direction have
been presented by Cromer [10], Ball [11, 12], and Kim
[13]. The results of Cromer and Ball also indicate that
suddenly applied axisymmetric pressures may couple breath-
ing mode response into preferred vibratory modes to cause
dynamic instabilities at pressure loadings lower than
static critical loads. These studies indicate that cou-
pling between parametrically excited flexural modes is
very weak and often may be neglected without compromising
the results.
Stricklin [14] and Wu [15] have made recent and
significant contributions to the nonlinear dynamic analy-
sis of shells by employing the finite element method. Re-
sults published by both of these investigators indicate
that accuracy comparable to methods emanating from the
equations of motion is obtainable if sufficiently complex
elements are used in large quantity. However, due to
longer computing times required for finite element tech-
niques, their major advantage clearly lies in the treat-
ment of complex structures where equations of motion are
not reasonably derivable. A seldom mentioned, yet almost
universal, appealing characteristic of finite element
techniques warrants a brief mention. This is that the
unknown dependent variables in a finite element solution
have a direct physical interpretation at all times as
the solution proceeds. Thus, to some extent, the analyst
never loses physical contact with his problem. This is
in contrast to most other analysis techniques applicable
to shell dynamics where the solution is carried forward
in some transformed mathematical space.
Klosner [16, 17] in a series of studies covering
a five-year period studied the effect of suddenly applied
axial loads on cylindrical shells. His solutions were ob-
tained by elimination of the circumferential coordinate
from nonlinear shallow shell equations and applying the
Ritz method or Hamilton's principle to the relations re-
sulting after introducing assumed solutions. His results
indicate that suddenly applied axial loads do not exhibit
the dynamic load factor effect found by others [9, 10, 11,
12] when suddenly applied normal pressures are considered.
These investigations reveal that a valuable indication of
the onset of dynamic instability can be found in the re-
versal in the trend of the time to first maximum deflec-
tion of the critical mode.
Another subcategory of solution methods in studying
dynamic shell stability by direct temporal integration is
known as the energy method. This technique circumvents
the derivation of the partial differential equations of
motion. The method involves deriving internal energy
expressions in conjunction with appropriate kinematic
relations, substituting assumed solutions into these
relations, and then applying Lagrange's equations to
the result. This procedure leads to either a coupled
system of nonlinear ordinary differential equations if
the spatial integrals are evaluated analytically or a
coupled system of nonlinear integro-differential equa-
tions if the spatial integrals are to be evaluated
numerically. In either case, the differentials are
temporal and must be solved numerically. Goodier [18],
Mclvor [19], and Lindberg [20] have presented solutions
to a particular type of dynamic instability of cylin-
drical shells utilizing this solution procedure. Hubka
[21] has presented alternate methods for numerical inte-
gration of the final equation sets for this problem.
This problem is concerned with the circumferential modal
coupling arising from the nonlinear membrane stiffness
terms in the equation sets obtained by this energy method
technique. The main results of these studies are that non-
linear coupling between circumferential modes may under
conditions of minor excitation or initial imperfections
in the preferred vibratory modes cause energy extraction
from the driven breathing mode into one or more preferred
modes such that the preferred mode becomes dynamically
unstable.
Utilizing this same energy method derivation,
Mente [22] has studied the dynamic nonlinear response
of cylindrical shells subjected to asymmetric pressure
loading. His computerized solution, known as DEPICS,
is able to handle a variety of coupled modes during the
temporal integration. The nonlinear stiffness is handled
by solving a linear eigenvalue-eigenvector problem at
each time step then adjusting the succeeding stiffness
value based on the nonlinear terms. Unfortunately this
rather extensive computer program has to date been unable
to satisfactorily define any regions of instability.
Significant Russian contributions to the study of
dynamic stability of shells have been made. Among these
are works by Shumik [23, 24] who also employed solutions
by direct temporal integration of expressions derived by
the energy method. These works deal with suddenly applied
loads using a linearized uncoupled dynamic stability theory.
His results indicate that cones with large taper buckle
dynamically at circumferential wave numbers higher than
the critical wave number for static loading.
Shiau [25] utilizing an energy solution studied
the effects of imperfections on the dynamic stability of
conical shells subject to axial impact. His results
agree with Klosner's [16, 17] in so much that a reversal
in the trend of the time to first maximum deflection of
the critical mode may be viewed as an indication of the
onset of dynamic buckling. The results, however, are
so critically dependent upon the choice of the single term
assumed buckling mode that the procedure appears to be
of limited general value.
A final division of solution methods in the cate-
gory of direct temporal integration is the method of
weighted residuals. The details of the procedure will
be covered in a later chapter. Nash [26] utilizing
shallow shell theory in conjunction with the Galerkin
method of weighted residuals investigated the response
of thin conical shells to dynamically applied axial
forces. He found that a rapid increase in end shorten-
ing may be used to denote the onset of dynamic instability.
Fisher [27] employed a similar solution method to
study the dynamic buckling of reinforced circular cylin-
ders subjected to suddenly applied axial compression load-
ings. He also found that a reversal of the trend of the
time to maximum deflection corresponds to onset of dynamic
buckling. Lakshmikantham [28] using shallow shell theory
in conjunction with smeared analysis and the Galerken tech-
nique investigated the same problem as Fisher [27]. His
solutions indicated that a jump in radial deflection was
the correct criteria for indicating dynamic instability.
An earlier work by Dietz [29] studying the same problem
by shallow shell theory and the Galerkin method indi-
cated that loss of stability is best defined as a sudden
drop in the axial load required to maintain a constant
end shortening rate.
By taking the results of the preceding literature
survey in context, several comments may be made concerning
the current state-of-the art in analytical studies of
dynamic shell stability. First, it appears that nonlinear
response problems which retain full circumferential modal
coupling yield the most reliable indicators of the onset
of dynamic buckling. Also the results appear to be highly
sensitive to theoretical approximations and to the quality
of the numerical solution procedures. Secondly, the study
of passage through bifurcation points is well developed
only for linearized dynamic stability techniques. The
few articles proposing to deal with nonlinear dynamic
stability unfortunately discarded the quadratic nonline-
arities, which are the most significant ones, retaining in-
stead only the cubic nonlinearities. This was done in order
to circumvent having to deal with inter-modal coupling. There-
fore, the results of these studies appear inconsistent.
It is the purpose of this dissertation to invest-
igate the dynamic stability of elastic shells of revolu-
tion. Two specific areas of this broad field are treated
in detail. First, this analytical study generates previ-
ously unavailable interaction relations for combined
dynamic loadings as they interact to cause passage from
a dynamically stable state to other states. Secondly, the
concept of linearized dynamic stability is extended to in-
clude geometrically nonlinear effects. One role of these
effects is to allow the possibility of autoparametric
excitation of preferred flexural modes by the driven modes.
The question of whether such excitation can occur during
the primary dynamic instability motion sufficiently to af-
fect the magnitude of critical dynamic loads is studied.
The theory derived herein is suitable for the study
of the dynamic stability of general shells of revolution.
In order to obtain numerical results, however, a specific
shape must be specified and the one used for numerical
analysis in this study is a truncated circular conical
shell structure. Conical frustum shells are in frequent
use as structural elements. Moreover, a conical shell
reduces to a cylindircal shell when the semivertex angle
of the cone becomes zero. Similarly, it reduces to a flat
circular plate when the semivertex angle of the cone ap-
proaches a right angle. Thus, the analysis and computer
program developed for cones can be generalized to handle
these two related problems.
A modified version of the subdomain method in
combination with circumferential modal analysis is de-
veloped for the solution method. A computer program is
constructed to obtain numerical results and the dynamic
stability characteristics of a conical frustum are stud-
ied under a variety of combined dynamic loadings. Inter-
action curves for static stability, linearized dynamic
stability, and nonlinear dynamic stability are generated.
A dynamic buckling criterion based on meridional mode
shape changes is established and the quantitative impor-
tance of nonlinear coupling effects is illustrated.
CHAPTER 2
THEORETICAL DEVELOPMENT OF NONLINEAR SHELL EQUATIONS
2.1 Shell Coordinates
A set of normal curvilinear coordinates (x, 0, z)
is chosen in the shell space such that a and B lie on the
undeforme d middle surface, z is perpendicular to this
middle surface, and z=O lies at the middle surface. The
equation of the undeformed middle surface is given in
terms of two independent parameters a and 0 by the radius
vector
r = r(a,B). (2.1)
To describe the location of an arbitrary point in the
space occupied by a shell, the position vector is defined
as
R(a,B,z) = r(a,B) + zin(a,B) (2.2)
where z measures the distance of the point from the corres-
ponding point on the middle surface along i The unit
vector i is normal to the surface.
2.2 Kinematic Relations
Acceptable nonlinear strain-displacement relations
for shells may be obtained from one of two basically dif-
fering approaches. One approach is to utilize the finite
strain tensor of the three-dimensional theory of elasticity
in conjunction with the Gauss-Codazzi relations for a sur-
face [30, 31, 32, 33, 34]. Another method is the two-
dimensional approach to shell theory that evades the ques-
tions of the approximations involved in the descent from
three dimensions. This approach defines strain as one-
half of the difference between the material and spatial
metric surface tensors. Sanders [35] presented the first
consistent nonlinear kinematic relations derived by this
method using a continuum called a Cosserat surface.
Naghdi [36] discusses both of the above methods in detail.
Nonlinear kinematic relations are seldom utilized
in shell analysis in their complete form due to their
complexity. Therefore, simplifications are routinely made.
Throughout the following analysis it is assumed that the
strains in the middle surface and that the rotations about
coordinate axes are small. These assumptions imply that
an element of area on the deformed middle surface is iden-
tical in size to an element of area on the undeformed
middle surface and that the difference between the
Christoffel symbols of the deformed and undeformed coordi-
nate systems is zero. These assumptions result in a con-
siderable simplification in the resulting equations with-
out unduly restricting the applicability of the solutions.
The strain-displacement equations shown below were
derived using the finite strain description of the three-
dimensional theory of elasticity in a manner similar to
Ogibalov [32] and the nonlinear terms were regrouped into
surface rotation expressions. The resulting kinematic
relations are
ea = (1/(l+z/R ))[U, /A+A, V/AB+W/R +(l/(2(1+z/R )))
x(-W, /A+U/Ra 2
e = (1/(I+z/R ))[V, /B+B, U/AB+W/R +(l/(2(1+z/R )))
x(-W, /B+V/R )2] (2.3)
e = (1/(2(1+z/R )))(V, /A-A, U/AB)
+(l/(2(1+z/R )))(U, /B-B, V/AB)
+(1/(2(1+z/R a)(1+z/R )))(-W, a/A+U/Ra)(-WB/B+V/R ).
In equations (2.3) the covariant components of the normal
and shear strains are denoted by e e e a respectively.
The components U, V, and W represent displacements of
an arbitrary point in the shell. R and R are the
principal radii of curvature and A and B are the coeffi-
cients of the first fundamental form. Comma denotes
partial differentiation.
2.3 Constitutive Relations
The relations between the components of strain
and stress in an orthotropic, linearly elastic material
are given by Hooke's law [37] as
e = oa/E -v a /E V a n /E +a T
a a / -vcn n c
e = -v a /E +cA/E -v n/E +a T
e. / n n /(
en = -v n oa/E -v n /E +,n/E +aT
(2.4)
2eB oa/G
2e = an/G
an an.
Young's modulus in the a and 3 directions is de-
noted by Ea and E respectively. The Poisson Ratio v(B
designates the contraction in the a direction caused by
a positive normal stress a in the direction. G is
the shear modulus in a plane which is tangent to the
(a,6) coordinate surface and (a(, a ) are the coeffi-
cients of thermal expansion in the a and directions
respectively. These relationships are valid for materials
subjected to stresses below the proportional limit. By
introducing Kirchhoff's hypothesis and allowing for the
symmetry of the elastic parameter matrix [38] the above
equations may be written as
aa = [1/(1-v O )][E e +v e -(a E +a E )T]
aB = [1/(1- a V )][E +v Ee -(a E +a E )T] (2.5)
_a = 2G e
2.4 Kirchhoff's Hypothesis
The reduction of the three-dimensional problem to
a two-dimensional one requires an assumption concerning
the variation of strain, displacement, or stress across
the shell thickness. To satisfy this requirement the
fourth assumption in Love's first approximation [39],
known as the Kirchhoff hypothesis, is introduced. This
hypothesis entails the vanishing of transverse shearing and
normal strains [33] and may be formulated as follows:
U(a,,z) = u(a,8)+z6a(a,6)
V(a,6,z) = v(a,6)+ze (a,()
W(a,B,z) = w(a,) (2.6)
e = u/R -w, /A 6 = v/R -w, /B
where u, v, and w are the components of displacement at
the middle surface in the a,4, and normal directions,
respectively, and 0 and 6 are the rotations of the
normal to the middle surface during deformation about
the B & a axes, respectively. The acceptance of this
assumption is due to its great clarity [40]. Although
Kirchhoff's hypothesis is a first approximation, its
applicability to nonlinear shell theory is well known
[33, 35]. The problem of determining the error intro-
duced by the hypothesis on the preservation of the normal
has thus far not been solved exactly [41]. Novozhilov
[42], Mushtari [43], and Koiter [44] have estimated that
the errors introduced by Kirchhoff's hypothesis are of
the order of h/R, which is small for thin shells.
Kirchhoff's hypothesis, modified to account for
transverse normal strain, has been employed in [31, 33,
34, 36, 45]. Effects of transverse normal strain were
not included in the present analysis.
2.5 Variational Principle
2.5.1 Hamilton's Principle
The Irish mathematician and physicist, Sir William
Hamilton (1805-1865), formulated his celebrated principle
in dynamics in which the governing equation depends ex-
plicitly on the energy of the system [46]. Hamilton's
principle is stated as an integral equation in which the
energy is integrated over an interval in time. In the
language of the calculus of variations, Hamilton's prin-
ciple states that the first variation of the time integral
of the difference between the kinetic energy T and the
potential energy P of a dynamical system is zero, that is,
t2 %
6 f (T P)dt = 0 (2.7)
t1
The equation is assumed to hold for all dynamical systems
whether they are conservative or nonconservative. The
quantity (T-P) is called the Lagrangian function and is
denoted by L. With this designation equation (2.7) be-
comes
t2
6 f Ldt = 0 (2.8)
t1
1
and Hamilton's principle then asserts that the first vari-
ation of the time-integral of the Lagrangian function is
zero. Hamilton's principle is employed in the next sec-
tion to derive the equations of motion and the admissible
boundary conditions for shells.
2.5.2 Equations of Motion
Of the many procedures available for the derivation
of the equations of motion for a differential shell element,
Hamilton's principle is superior in nonlinear shell theory
due to its efficient treatment of problems involving curv-
ilinear coordinates and because it gives the admissible
boundary conditions, natural or forced, that are to be used
with the theory [47, 48]. With the total potential energy
P expressed as the difference between the internal strain
energy and the potential energy of the external forces, we
obtain from equation (2.8)
t
6f {f [(po/2)V i-(1/2)siJij]dv+f si Vdsldt = 0 (2.9)
t 0 0s
Here the density of the undeformed body and the Lagrangian
components of the displacement vector are p and V., re-
spectively. The symmetric Cauchy stress tensor sij is
measured per unit area of the undeformed body and must
be referred to the base vectors in the deformed body and
the strain tensor, Yij, is defined by equation (2.3).
The components of the surface force ss are referred to
the base vectors in the undeformed body and the volume
and surface integrals in equation (2.9) must be extended
over the volume and surface of the body in its undeformed
state.
Introducing the kinematic, constitutive, and
Kirchhoff relations from previous sections, we eliminate
the time derivatives of the variations by integrating by
parts with respect to time and require the virtual dis-
placements to vanish at the end points of the arbitrary
interval t1 t t2, thus obtaining the equations of motion,
stress resultants and couples, and boundary conditions.
The equations of motion derived in this section
reflect the assumptions of Kirchhoff's hypothesis. Small
strains and moderately small rotations were also assumed.
Thinness assumptions were delayed until after the equa-
tions were synthesized.
The equations of motion are:
(BNO), +(AN(), +A, N B-B, N +ABQc/R -AB(NO +Na O )/Ra
+ABpa = ABphu(l+h2/12R R )+ABph3(1/R +1/R ) /12
(BNaB), +(ANB), +B, Naa-A, Na+ABQ /R -CAB(N 3 +N Be )/R
+ABp1 = ABphV(I+h2/12R R )+ABph3(1/R +1/RB1 6/12
(BQa), +(AQ3), -AB(Na/R +NB/R )-(Be N+BN (2.10)
-(AG N a+AO^ N), +ABpn = ABphW(l+h2/12R R )
(BMO'),+(AMB1), +A, MaB-B, M -ABQa
= ABph3[(l/R +l/R 3 )U+]/12
(BMW3), +(AM3), +B, MBa-A,1Mc-ABQ3
= ABph3[(I/R +I/R )V+ ]/12
where (") denotes second partial differentiation with
respect to time and p p, pn represent the applied
surface loads.
The stress resultants and stress couples are given by
Nc h/2 a
{ = f { 1(1+z/R )dz
N a -h/2 ao0
{N }
N
{Mm }
M m
M
M }
QO
N
HS
h/2 B
= f { 1(1+z/R )dz
-h/2 az U
h/2 a
= f h{ }(1+z/R )zdz
-h/2 0c
h/2 a
-h/2 a}(1+z/R )zdz
h/2 n"
= / {0 1+z/RB}dz
-h/2 o5n 1+z/R
(2.11)
h/2
= f o~dz
-h/2
h/2
= f a zdz
-h/2
where h is the shell thickness.
Figure 2 illustrates the relations between the
stress resultants and couples to the generalized coordi-
nates and to the shell element.
Z,W
M B
C, u
Figure 2. Shell Element
7
2.5.3 Boundary Conditions
From equations (2.9), (2.10), and (2.11) we get
the following natural boundary conditions that arise
from the requirement of force balance. Along the edge
of constant a the boundary conditions are
N = N or u = u
N +M /R = N or v = v
(2.12)
Q U+Ma-N -N6 = Qa or w = w
Ma = Ma or 6 = 6
The boundary conditions along constant B may be obtained
from equation (2.12) by interchanging a with B and u with
V.
2.6 Synthesis of Equations
The usual procedure followed is to reduce the
number of equations and unknowns to a more manageable
number by eliminating Qa and Q from the equations of
motion (2.10) which are thus reduced from five to three.
The force and moment resultant expressions (2.11) to-
gether with the constitutive relations (2.5) are then
substituted into the equations of motion. Finally, the
strain-displacement equations (2.3) and (2.6) are substi-
tuted, to yield three differential equations of motion
having u, v, and w as dependent variables and a, B, and
t as independent variables.
2.6.1 Shells of Revolution
In the preceding analysis Ra and R are arbitrary,
but in all subsequent applications these principal radii
will refer only to shells of revolution. It can be shown
[47] that for a shell of revolution the principal radii
may be expressed as
Ra = -[l+(DR /x3) 2 3/2 /2R /x3 2)
(2.13)
R = R /l+(DR /Dx3)2
where R is the radius of the latitude circle and x is
the coordinate coincident with the axis of revolution.
2.6.2 Conical Shells
The above described synthesis may be carried out
using equation (2.13); however, the equations would be
extremely unwieldy. Thus, the procedure will be followed
only for a specific shell of revolution shape and the one
used in the remainder of this study is a circular conical
shell. A conical shell reduces to a cylindrical shell
when the semivertex angle becomes zero. Similarly, it
reduces to a flat circular plate when the semivertex
angle approaches a right angle. Therefore, the following
analysis can handle cones as well as the two related
problems. The conical coordinate system adopted is shown
in Figure 3. The coordinates, coefficients of the first
fundamental form, and principal radii become
a = s 8 = 0
A = 1 B = s sin a (2.14)
R = c R = s tan a.
In all subsequent relations, a will denote the semivertex
angle.
Assumptions relating to displacement magnitudes are now
introduced into the analysis. The squares of inplane
displacements are assumed to be small with respect to the
squares of normal displacements and the shell thickness
is assumed to be small compared to the radius.
A
f--
A -A
Figure 3. Conical Shell Coordinates
Neglecting rotary inertia terms, the nonlinear
equations of motion are
Ns+sNs -N +N,0 = sphu
-ctn2 vN /s+ctna w, N /s+N, +2Ns+sNs+ctna w,sNs
0 so sO
+ctna M, /s+2ctnc M /s+ctna Ms = sphV
w N5s N0 0 w
w, sN+sw, ssN +sw,sNs-ctna N -ctna v, N /s+w,W N /s
0 0 sO se sO
-ctna vN, /s+w, N, /s-ctna v,NSN +2w, sN -ctna vN's
sO sO s s 0 0 s0 sO
+w, Ns +w, N, S+2Ms +SMss-Ms+M, /s+2M, /s+2Mse +sp
5'i 's s 5 s Iss- s p / s' P
= sphW
where
(2.15)
= 0 sina.
The nonlinear stress resultants and stress couples
are
Ns = [Eh/s(l-v 2)]{su,s+s(w,)2/2+v[v, +u+ctna w+(w, )2/2s
-ctna vw, /s]}
Nso = Nas = (Gh/s)(sv, +u, -v+w, W, -ctna vw, )
N6 = [Eh/s2(1-v2)]{sv, +su+sctna w+(w, )2/2-ctna vw,
+ctn2 v2/2+v[s2us+s2(w,s)2/2]} (2.16)
Ms = [Eh3/12s(l-v2 )]{-sw,ss+ctna u, s+ctna(w,) 2/2
+v[-w, 0/s-W,s+ctna(w, )2/2s2]
Mso = M -s = (Gh3/12s)(-2w, s+2w, /s+ctna w, sw, /s)
M0 = [Eh3/12s2(1-v2)](-w, -sw, -ctna u-ctn2 a w-vs2w,) .
Substituting equations (2.16) into equations
(2.15) yield three equations of motion. Noting that
Gh = K(I-v)/2
we have
K{[-u/s+u,s+su, ss+u, /2s-3v, /2s+v, s/2-ctna w/s
+(w,s/2s)(w, -ctna v, )+(w,s/2s)(w, -ctna v)+(w, )2/2
-(w, )2/2s2 +ctna vw, /s2-ctn2a v2/2s2+swsw,ss
(2.17)
+v[-u, /2s+v, /2s+v,s/2+ctn w,s-w, s/2(w, s+w,/s
-ctna v, /s)+w, /s(w, s/2-w, /2s-ctna v,s+ctna v/s)
-ctna vw, s/2s]} = phsu
K{3u, /2s+u, s/2-v/2s+v,s/2+sv, ss/2+v, /s+ctna w, /s
+w, (w, ss/2+w, /s2)+W, (w,+w, /s)/2
+ctna[-sv(w,ss/2+w, /s2)+W,s(-v+u, /2)+w, u/s]/s
(2.18)
+v[-u, /2s+u, s/2+v/2s-v,s/2-sv, ss/2+w,sW, s/2-w,sw, /2s
-w, ssw,/2+ctna u, s/s(-ctna v+w, )+ctna w,s(v-u, /2)/s
+ctn vw, ss/2]} = phsV
K{-ctna (u+v, +ctna w)/s+ctna[(w, )2/2+ww, ]/s2+su,s W,ss
+Ws(Uss+Vs/ '^ ^^ 't 9 9N ss
+w, (su, +v, /2+u,w /2s)+w, v, /s 2+w, (v, 4s 2+V,ss/2
+u, '/2s)+w,'s(u, /s+Vs)+u,s'W, +uw, /2+u 2 w,-/2s2
-vw, s/s-v, sW, /2s+vw, /2s 2-v, W,s/2s+ctnct(-v,2 /s2-uv, /S2
-vv, /s2-u, v/2 's -V,/2-uV, 's/2s+vv, s/s-v, ss/2-u, sv/2s
-v 2/2s2 )+ctn a2 (-v, w/s 2)+v[-ctna u, +w, s(v, +ctna w)
+w, (U,S /2s-v, ss/2)+us W, 4/s+w, s(Vs/2+ctn ws/2-u, w/2s)
+uw,+ss+U,sW, s-V, sW,s-u,W,'s/s+vw, /s+u, w, /2s2 (2.19)
+v,s W,/2s-vw, /2s2+v, Ws/2s+ctna (-u, sV,/s-us,sv/2s
+v,2/2+u, v, /2s+vv, /2-u ,v/2s2-vv, /s+v2/2s2)]}
+D{-swssss-2wss2 /s-w, /s3+2w'sq/s2-2w'sss-4w, /s3
'ss s s s^s ss s /sss ss, -4
+W,'s W, /s2-'sW, /s3+Wss W, /s2+(w, ) 2/s2
+v[(ctna/s2)((-7/2s)(w, w,'s)+3(w, )2/s2+ww,s W,/2s
-WssW,-W,'s W,'s)]}+spn = phs
where the membrane stiffness is
K = Eh/(1-v2)
and the flexural stiffness is
D = Eh3/[12(1-v2)1.
(2.20)
(2.21)
CHAPTER 3
NONLINEAR DYNAMIC STABILITY
EQUATIONS OF CONICAL SHELLS
The nonlinear dynamic stability equations are
obtained from the nonlinear equations of motion by
employing a perturbation analysis. Substituting the
displacements
u = uI + uB
v = vI + vB
w = I + wB
(3.1)
and the stress
Ns N s + Ns
= NI B
N = N + N
N N + N
I B
resultants and stress couple expressions
M1 = MS + Ms
I B
M = M1 + MB
Mso = Ms + Mf
I B
(3.2)
and the load relation
P = P + PB (3.3)
into equations (2.15) results in two coupled sets of equa-
tions. The quantities with subscript I are measured from
the undeformed state and the quantities with subscript B
are small but finite perturbations. This gives one equa-
tion set describing the motion in the stable space surround-
ing the undeformed state and another set, the stability
equations, describing the motion during passage thru the
nearest bifurcation point.
The resulting set of nonlinear dynamic stability
equations are as follows where the subscript B has been
dropped for convenience:
K[(-u/s+u,s+su,ss+u, /2s-3v, /2s+v, s/2-ctnc w/s)
+v(-u, /2s+v, /2s+vS/2+ctna w,s)
(3.4)
+(ww/2s+w2 2 2)
+(WsW, /2s+w's, w, /2s+w,s/2-w, /2s +sw,sW,ss
(-w /-w /sww /2s-w2/2s2)] = phsU
Ws/ sW, 2 w w
K[(3u, /2s+u, s/2-v/2s+v,s/2+sv, ss/2+v, /s+ctna w, /s)
+v(-u, /2s+u, /2+v/2s-v,s/2-sv, ss/2)
(3.5)
+(w, W,ss/2+w, w, /s2+wsW, /2+w, sW,/2s+ctn2a ww, /s2
+v(w,s W, /2-w,s W,/2s-w,ssW, /2)] = phsV
K{-ctna(u+v, +ctna w)/s+,v(-ctna u,)
+ctna [w/2s 2+ww, /s2+V(w, /2+ww, )]}
+D{(-swssss- 2wssq /S-W, O /s3+2wsw /s2-2sss
-4w, /s3+W, ss/s-W, s/s2)+ctn[w, 2+w, ssw /s2
+ws 2/S2 +W 5sss W,'s5/s3-ws '4W,/s3+w, Wss/s2 (3.6)
+W,'s W, /s2+v(-3w, w's /s +3w, /2s +Ws W, /s3
-w,ssWw/s2-, s Ws/s2)]} + spn
+w,sNs+sw, ssN+sw, sNIs wN/s+w,l N /s = phsw .
The single underlined quantities represent non-
linear membrane stiffness terms. The double underlined
quantities represent nonlinear bending stiffness terms.
The nonlinear stability terms are denoted by a dashed
underline. These equations embody the assumptions that
wI < wB
and that
uB < w
vB < WB
In the present study, a thin truncated conical
shell is loaded by a constant axial load and a nearly
axisymmetric Heaviside pressure load. The membrane-
state stress resultants are given by
NS = -P/(Trs sin2a)-pn(stanoa)(1-cosxot)/2
NI = -pn(stana)(l-cosx t)
where P is the applied axial load, pnI is the magnitude
of the Heaviside pressure load, and A is the natural
frequency of the breathing mode of the shell.
CHAPTER 4
METHOD OF SOLUTION
4.1 Circumferential Modal Analysis
The nonlinear partial differential equations
derived in previous chapters are spatially two-dimen-
sional. The two directions are meridional and circum-
ferential. In this analysis, the two-dimensional equa-
tions are reduced to one-dimensional equations by
utilizing circumferential modal analysis. The inde-
pendent variable 0 is eliminated by expanding all
dependent variables into sine or cosine series in the
circumferential direction. As a result of the trig-
nometic series expansions, there is one set of govern-
ing equations for each circumferential wave number
considered. For a linearized analysis the sets are un-
coupled. For a nonlinear analysis the equation sets are
modally coupled through the quadratic terms. This cou-
pling arises when use is made of the trignometric prod-
uct identities
cos(ie)cos(je) = (1/2)[cos(i-j)6+cos(i+j)6]
(4.1)
sin(ie)sin(je) = (l/2)[cos(i-j)e-cos(i+j)e]
40
in order to facilitate the equating of like coefficients
of Fourier expansion terms which is required to identify
the equation sets.
Substituting the expressions
u(s,i,t)
v(s,p,t)
w(s,4,t)
p(s ,t)
Q
n=o
Q
= z,
n=o
Q
= z,
n=o
Q
= o
n=o
u (s,t) cos(np/sina)
v (s,t) sin(ni /sina)
Wn (s,t) cos(ni/sina)
p (s,t) cos(n4/sina)
into the nonlinear dynamic stability equations (3.4 -
3.6) and equating like coefficients of Fourier expan-
sion terms yield the following:
(-l/s){l+[n2(l-v)/(2sin2 )] u +u n's +su n'ss
+(n/2s sina)(-3+v)v +(n/2sina)(l+v)v n -(ctna w /s) (4.3)
+vctna win's = phsu n/K
(4.2)
(n/2s sina)(-3+v)u (n/2sina)(l+v)u n's
+{-[(1-v)/2]-(n2/sin2 )}(v /s)+[(1-v)vn's/2] (4.4)
+[(1-v)svn',ss/2]-[n ctna w /(s since ] = phsV n/K
-ctna [(un/s)+vun's+(n/sina)(v /s)+(ctna wn/s)]
+(D/K){[(4n2/sin2a)-(n4/sin4 )](w /s3)+[(2n2/sin2a)+1]
x[(Wn'ss/s)-(Wn's/S2)]-2wn'sss-swn'ssss}+ p n/K
-[P w /(KTrsin2a)]+(1-cosX t)[-s tana B ,s
I n 'ss /2)+(2n c2 s
-(s2tana n, /2)+(2n2 n/sin2a)-(2n n/sina)]/K
c s cs/4) (4 5)
+nc{ctna{(-n2w2/4s2sin2a)+v\)[(w ,2/2)+(w n,2/4) (4.5)
+w w +(w w /2)]}+(h2ctna/12){w 2 +(w 2/2)
+o o'ss n n'ss O'ss n SS
-[n2/(2s2sin 2a)](nw n ss-w ,2)}}+n {ctna[(-n2 wo /
(s2sin2a))
-(no's n'sW oss n o n'ss 'ss nss
-(n2/s2sin2a)wo'ss Wn]} = phsOn/K
where
n = circumferential wave number
S= 1 n=o
S o n>o
n o
n zw2c c c o n=o (4.6)
c1= ( Pn) i-n P I n>o
1= w./2[ns p( i+n)+ sli-n p i nl],
s i=-0
s n=o
n>o
s = l in+l
p =o i=n
-1 isn-1
4.2 Elimination of Meridional Derivatives
The next phase of the analysis involves the re-
duction of the partial differential equations sets result-
ing from the circumferential modal analysis to ordinary
temporal differential equation sets which are then amen-
able to numerical solution.
Referring again to Figure 1 it is seen that the
available solution methods for this task fall generally
into four categories. These are finite element methods,
finite difference methods, variational principles, and
methods of weighted residuals. As noted by Finlayson
[49], for certain types of linear problems these methods
can be shown to be equivalent to each other. From the
literature abounding on shell dynamics one can safely
conclude that any of these methods when judiciously ap-
plied will give acceptable results. There are, however,
marked differences between the methods in the amount and
quality of manual mathematics and the quantity of machine
computing time required to give equivalent results.
The desired characteristics of the solution method
for the particular problem of interest in this disserta-
tion are now discussed. The solution method utilized
should yield an approximate direct solution where con-
vergence would be monotonic and assured. Spatial control
of error residuals was desired. Because the taper in the
conical shell results in nonsymmetric stiffness matrices,
a solution technique applicable to non-self-adjoint prob-
lems was required. Finally, it was hoped that the chosen
solution method would result in short computer compile
times because experience has shown that the majority of
computer time required for such a problem is in successive
compile times during development and debugging.
In matching these desired characteristics to the
known characteristics of the classical solution methods
an elimination process was begun which finally resulted
in the algorithm used herein.
Finite element techniques were eliminated from
consideration because this method is not competitive with
other methods in problems with regular geometries and
formulatable boundary conditions lying along coordinate
lines. This method requires significantly more manual
mathematics and computer execution time than equilibrium
equation solution methods in order to give comparable
results.
Finite difference methods have been shown to be
extremely powerful in solving nonlinear dynamic stability
problems involving shells. This method has been used with
particular success when coupled with circumferential mode
analysis. It was shown, however, by Radosta [50] that
convergence is slow when a meridional mesh is used with
a shell under high axial load. In the particular problem
at hand, i.e. a conical frustum subject to a high axial
load combined with a time dependent pressure load, it was
felt that a suitable selection of number and spacing of
meridional mesh points would degrade into a highly time
consuming trial and error effort.
The final phase of solution method selection con-
sisted of choosing one of the four methods of weighted
residuals: Galerkin, Least Squares, Collocation, and Sub-
domain. The least squares method was not seriously con-
sidered because it yields ordinary differential equations
of the second degree in the temporal operator. The
collocation technique was eliminated because it enhances
the same point location selection drawbacks as does the
finite difference method. Again, this is normally an
insignificant problem but was made significant here by
the expected buckled states resulting from the combined
type loadings being considered. The final choice between
the Galerkin and subdomain methods was based on the method
of error residual control. In the Galerkin technique the
error produced by each individual term in the assumed solu-
tion is forced to average to zero over the entire length
of the shell. In the subdomain method the error in the
total assumed solution is forced to average to zero over
each subdomain. This type of spatial control for error
residual was considered superior to the modal type of
error control offered by the Galerkin method and there-
fore, the subdomain technique was selected for use in this
investigation.
4.2.1 Subdomain Method
A description of the theoretical basis of the sub-
domain method is now in order. Consider the equation
L[u(x,t)] = 0 in V (4.7)
where L is a nonlinear, non-self-adjoint, ordinary
differential operator. If the approximation solution is
represented in the form
m
= Z ai (t)4i(x) (4.8)
i=l 1 1
and this is substituted into equation (4.7), we obtain
an approximation
L[u(x,t)] = e (4.9)
where E is the error. We now impose constraints on E by
requiring it to average to zero over a subdomain v of
the total domain V
fL[u(x,t)]dv = 0 (4.10)
v
The number of subdomains v are chosen to equal m the
number of time dependent unknowns in (4.8). The differen-
tial equation, integrated over the subdomain is then zero,
hence the name subdomain method. As m increases, the dif-
ferential equation is satisfied on the average in smaller
and smaller subdomains, and presumably approaches zero
everywhere. This process yields a dynamically coupled
second order set of nonlinear temporal differential equa-
tions which may in turn be solved by numerical methods.
Equation (4.8) is called the trial function and
may be selected in a variety of ways. It may be chosen
to satisfy the boundary conditions, but not the differen-
tial equations or to satisfy the differential equations
but not the boundary conditions. In order to simplify
the required hand calculations, a set of complete power
series were chosen as the trial functions for this prob-
lem. The trial functions were then forced to satisfy
the boundary conditions by a modified application of the
Lagrangian multiplier method which can be more properly
thought of as a coordinate transformation to require an
assumed power series solution to satisfy the essential
boundary conditions. An approximate satisfaction of the
differential equations was then obtained by the subdomain
method.
4.2.2 Method of Satisfying Boundary Conditions
A description of the technique utilized to satisfy
boundary conditions will now be presented. Consider the
trial function (4.8) in matrix form,
= L J {a} (4.11)
(1xm)(mx1)
and the j essential boundary constraints,
[C ]{a} = 0. (4.12)
(jxm) (mx1)
Equation (4.12) may be partitioned as follows
[ C1 C2 ] 1 q 1
[jx(m-j)] (jxj) q2
(mxl)
Equation (4.13) may be rearranged as
{q2} = -[C2 -I [ C1 ] {ql}
(jxl) (jxj) [jx(m-j)][(m-j)xl]
and noting that
{ql} = [I]{q }
we obtain
{a} = [ B ] {ql}
(mxl) [mx(m-j)][(m-j)xl ]
and (4.11) becomes
u = LpJ[B]{qI}
where [B] may be considered as a projection operator.
(4.13)
(4.14)
(4.15)
(4.16)
(4.17)
49
After nondimensionalizing equations (4.3 4.5)
according to
h = s2h
S = P P
I cr I
s = s2s
S= s2u
Wn = 2n
-Sw
P
cr
= 2rEh2cos2a//3(l-v2)
(4.18)
T = /E/p(1-v2) t/s2
Pn = Eh/(l-v )s2 Pn
\ = s2/p(1l-2)/E An
derivatives are likewise
nondimensionalized.
the trial functions are introduced into the equations.
The trial functions are
un an.
"n ni
v N- b .
n = ni s (4.19)
w io Ci
Pn Pni
After applying the transforms satisfying the
boundary conditions, the equations are then integrated
over (N-j) subdomains and the average residual error for
each subdomain equated to zero. During the original
formulation it was noted that this approach leads to a
singular mass matrix for equally spaced subdomains. To
circumvent this problem the equations were mutliplied
through by s. This is similar to a procedure known as
the method of moments first employed by Polhausen [51].
Having completed the above sequence of operations
results in Q sets of nonlinear second order temporal dif-
ferential equations. Each set contains 3 (N-j) equations
and the sets are coupled in the nonlinear terms. In ma-
trix form a typical set may be expressed as
[m][B]{q}(n)+[k][B]{q}(n)+[D][B]{q}(n)+{Nl.}(n) {p}(n)
(4.20)
where
[m] mass matrix
[B] matrix requiring power series solution to
satisfy the prescribed boundary conditions
{q}(n) generalized coordinates of the nth cir-
cumferential mode
[k] stiffness matrix
[D] stability matrix
NL. Lq J[B ]T[C ][B ]{q }- power series products
{P} load vector .
51
As an example, for a problem retaining two circumferen-
tial modes with eight subdomains along the meridian,
equations (4.20) would represent a total set of forty-
eight coupled nonlinear second order ordinary differen-
tial equations to solve.
CHAPTER 5
NUMERICAL ANALYSIS FOR CONICAL SHELLS
The sets of equations (4.20) were coded into
computer language for numerical solution. A particular
shell specimen was selected for study. The salient
geometry and boundary conditions for the shell studied
are given in Appendix I. Fortran IV, level H machine
language was used and the program was constructed for
double precision arithmetic. A flow diagram of the oper-
ational characteristics is given in Appendix III. A list-
ing of the program is given in Appendix IV.
5.1 Functional Description of the Computer Program
If inplane inertia is considered ignorable with
respect to normal inertia, then equations (4.20) separate
into coupled sets of 2(N-j) nonlinear algebraic equations
and 2(N-j) nonlinear first order differential equations.
The flow of computations proceeds as follows. The
linear and nonlinear stiffness matrices, mass matrices,
boundary transform matrices, buckling matrices, and load-
ing vectors are generated by the program. At time equal
to zero the problem solution begins with known external
loading and initial conditions. The known normal deflec-
tions are substituted into the first two equations of
equilibrium and the resulting u, v vectors are generated
utilizing a Gauss elimination technique with complete
pivoting. If nonlinear effects in the first two equa-
tions of equilibrium are considered, this mode is itera-
tive in nature. Having generated the u's and v's for
time zero, they are substituted into the sets of 3rd
dynamic equilibrium equations for temporal integration.
The time integration is accomplished by a Hamming pre-
dictor-corrector integration technique. This generates
the w's for the next time step. Normal velocities are
also computed. This completes the sequence of advance
to the first time increment beyond time zero and provides
the required initial conditions in order to begin the
next advance in time. The sequence then recycles and the
temporal advance continues. At any time interval the an-
alyst chooses the solved for dependent variables are
transformed from the solution space to real space and
displayed as output. Each output consists of the com-
plete meridional profile for each mode in the solution.
The build up of errors in the numerical solution
is controlled by the following mechanisms. A loss of sig-
nificance control indicator is incorporated into the Gauss
elimination routine which warns the analyst when solutions
are being derived from a nearly singular coefficient
matrix. An error bound control is available in the
predictor-corrector integrator which limits this type
of error by automatically adjusting the size of the
time step. Round off errors are minimized by making
all calculations accurate to sixteen significant figures.
Such rigorous control of possible numerical errors is
essential in dynamic stability analysis because of the
possibility that numerical instability might be mistaken
for actual structural instability.
5.2 Checking the Computer Program
The control of clerical and conceptual errors in
the construction of the computer program was done by con-
tinually checking back to simpler known cases. In some
instances the known cases could be obtained by hand cal-
culations, in others, the results of other analysts were
used. For the shell specimen studied, the computer pro-
gram has been verified by checking against known solutions
for static stability, static deflection, free vibration
modes and frequencies [39], linear dynamic response, and
regions of parametric instability [52]. A convergence
study was also conducted to determine the convergence
properties of the solution for increasing numbers of
subdomains. The results of this study are presented in
Appendix II.
5.3 Dynamic Stability Results
After the above verifications were completed, the
program was used to study the linearized and the nonlinear
dynamic stability of a short conical frustum with simply
supported boundaries subject to various levels of constant
axial loads and critical levels of nearly axisymmetric
Heaviside pressure blasts. Eight subdomains along the
meridional direction were used. The driven breathing mode
and an arbitrary preferred flexural mode were retained in
the analysis. Complete coupling between the modes was re-
tained through the quadratic terms for the nonlinear case.
The critical static normal pressure for the speci-
men studied is 42.8 pounds per square inch and the critical
circumferential wave number is eight. This represents a
mid-span deflection in the breathing mode of 9 percent
of the thickness. The ordinate scale (W) of all figures
depicting time history responses is normalized by this
critical deflection. The critical static axial load is
40,400 pounds. The applied loads for time history graphs
will be noted as
(a, b, c)
56
where
a = P/Pcr, b = p /p cr, c = p /pcr.
Additionally, note that
Po = POH(t)
Pn = Pn H(t)
where
H(t)
denotes a Heaviside pressure pulse.
5.3.1 Dynamic Modal Preference
The program was used to determine the critical
flexural modes during dynamic instability. The results
are shown in Figure 4 and the critical static buckling
modes are indicated. Because only a few closely spaced
critical modes dominate the flexural response and because
I
coupling between these modes is very weak [10], the param-
eter study that follows is reduced by considering the in-
stability of only a sequence of single flexural modes to
find the mode indicating instability at the lowest load.
5.3.2 Linearized Dynamic Stability
The linearized dynamic analysis was obtained by
considering the vector {NL}(n) in equation (4.20) to be
zero. This means that the internal loads from the pre-
buckling state may interact with postbuckling deflections
to cause primary dynamic instability (i.e. snap through)
but internal parametric instabilities due to the pulsat-
ing nature of the membrane state are prohibited.
The critical circumferential wave number for
dynamic instability was found utilizing this analysis and
the results are shown in Figure 5.
While conducting a series of linearized dynamic
stability computer runs at various combinations of external
loadings, a feature of the numerical model was discovered
which represents an improvement over other computer programs
used for this purpose. This feature is that the meridional
profile undergoes a marked and rapid change of character
during runs at or above the critical dynamic load. For
loadings below critical, the meridional profile pulsates at
I
-0 -0
0 (1)
owo
E r-
4-) 3 (3
Co -=
one half wave, slightly skewed due to the conical taper.
At or above the critical load the meridional profile
begins to pulsate as a half wave but later snaps to two
or more half waves. The time of snap, number and shape
of the waves in the buckled meridional profile depend on
the axial load to normal pressure ratio and their rela-
tion to static critical loads. Examples of this are shown
in Figures 6 and 7. It is felt that this improved ability
of detecting dynamic instabilities results from two areas
of the analysis. One being the spatial error residual con-
trol offered by the modified subdomain method employed and
the other being the strict control of error bound in the
numerical work. Additionally, it should be noted that
most published works on this type of problem utilize one
term assumed solutions for preferred modes and clearly
could not see this effect in their analysis.
A large number of linearized dynamic stability
runs were made at various combinations of static axial
and Heaviside pressure loads. The results were collated
into a dynamic stability interaction curve and are shown
in Figure 10.
5.3.3 Nonlinear Dynamic Stability
In the case of nonlinear dynamic stability, the
possibility of autoparametric excitation of preferred
6]
-0
m
0
4-
4-)
a, CD V
E
u- *
E
(Nj -0
S.-
.:I- CNJ
o 9
00 C
C-
0
00 0
V) E -C -C
o o
o o
Co C
C) C
cO o
SS-
4-
o
-4.
0 ---
\ E
CMt-~
5-
-,,
MENNEW
flexural modes by the pulsating breathing mode is allowed,
in addition to snap through or primary dynamic instability.
The question here is whether such excitation will occur
during the primary instability motion of the preferred mode
sufficiently to affect the magnitude of the external loads
required to induce dynamic instability. In contrast to the
external parametric excitation problems, the total energy
of the shell for the present problem is constant for all
t>0. Thus, the unstable modes when growing must derive
their energy from the breathing mode with a corresponding
decrease of energy in that mode. To evaluate this inter-
play a full nonlinear analysis is necessary.
Conceptually, the third equation of nonlinear dy-
namic stability may be expressed in an average sense as
follows:
2 2
w + w + k(u0 ) + f(w) p H(t)
(5.1 a,b)
[Xn g(w )]wn + 1(n ,vn) + h(w )wn = p H(t)
where u w0 represent breathing mode deflections measured
from the undeformed state and u v w represent the
poststable growth of the buckling mode. X and n denote
the linear eigenvalues of the respective modes and f, g, h,
k, and 1 are differential operators. The superscript a
denotes an average over the shell length. Since w0 is
excited by a Heaviside pulse, the w0 response is initially
periodic and equation (5.1 b) takes the character of an
inhomogeneous Mathieu equation. The primary buckling
terms denoted by a single underline cause an insignificant
shift to the left on a parametric stability diagram. In
the absence of imperfections in the preferred flexural mode,
denoted by the double underlined term, equation (5.1 b)
lies well within a stable zone of a Mathieu diagram and no
significant parametric growth of the buckling mode would
be expected.
For finite deterministic imperfections, however,
the location of stability boundaries becomes more difficult
and for significantly large p H(t), solutions to equation
(5.1 b)may exhibit parametric growth at or below the dynamic
loads required to produce immediate snap through and may
lead to a delayed dynamic snap through. Whether or not
this beating type of growth occurs with sufficient quick-
ness to have an effect on the primary dynamic instability
can be determined only by direct temporal integration of
two coupled sets of equation (4.20).
A series of runs with the computer program was made
to determine the significance of this effect. Typical re-
sults are shown in Figures 8 and 9. As in the case of
t.0
cy')
H~
cy*) C:
C
0 0
o 0
0 0
o 0
C \J
o o
LO U)
0 0
OO)CY
COOc
moom
linearized dynamic stability, a large number of computer
runs were made at various combinations of static axial and
Heaviside pressure loads and the results collated into
a nonlinear dynamic stability interaction curve shown
in Figure 10.
The linearized dynamic stability curve lies below
the static stability curve in Figure 10 due to the in-
creased severity of a suddenly applied load over that of
a statically applied load. The peak membrane forces de-
veloped by the breathing mode which is excited by a Heavi-
side pulse are twice as large as that of an equivalent
static case and contribute more to the developemnt of
instability in preferred flexural modes than a static
situation. This is known as the dynamic load factor effect.
The nonlinear dynamic stability curve lies below the
linearized dynamic stability curve because internal vibra-
tions interact between coupled modes so as to produce un-
stable beating resonances in the preferred buckling mode.
Such unstable vibrations lead to a delayed dynamic snap
through and this phenomenon is called autoparametric ex-
citation.
Static Stability
Linearized Dynamic
.... Stability
-_.=_ Nonlinear Dynamic
Stability
n = 8
0 0.5
p H(t) / pc
cr
Dynamic Stability Interaction Curves
1.0
0.5
0
1.0
Figure 10.
CHAPTER 6
RESULTS AND CONCLUSIONS
In the present analysis, a nonlinear shell theory
is derived and employed to study the nature of nonlinear
dynamic instability of a truncated thin circular conical
shell structure which is considered to be loaded by a
constant axial load and a nearly axisymmetric Heaviside
pressure load. A solution method is developed which sat-
isfies a boundary condition exactly and which converges
toward the exact solution of the governing equations with
increasing subdomains. The method avoids the necessity
of assuming the shapes of prebuckling or postbuckling
meridional profiles.
The results of this study confirm the feasibility
of the method of solution developed for this analysis. The
method converges rapidly yet can be applied to nonlinear
problems with minimum amounts of manual mathematics. The
trial functions are simple to manipulate yet satisfy any
formulatable boundary condition. The method is particu-
larly suitable for problems where poststable modes are not
known in advance as is usually the case for combined dynamic
loadings.
For the particular loading conditions studied,
the results indicate that dynamic instability will occur
at loads below critical static loads. This reduction in
critical dynamic loads was shown to be the result of both
the dynamic load factor effect and of autoparametric ex-
citation. The interaction curves which are generated
illustrate these effects quantitatively under conditions
of combined dynamic loading.
A criterion for dynamic buckling is established
based on meridional mode shape changes. This ability to
detect sudden jumps in the meridional profile aids to
verify instability detected by divergence on displacement
time history curves and provides additional information
about the poststable state.
APPENDIX I
APPENDIX I
SHELL SPECIMEN AND BOUNDARY CONDITIONS
The conical frustum utilized in this investigation
has the following properties (see Figure 3):
Material -
a = 200
sI = 5.85
s2 = 14.22
R1 = 2.0
R = 4.863
h = 0.02
v = 0.3 .
1020 steel
inches
inches
inches
inches
inches
The boundary conditions assumed for the analysis
are:
u = v = w = M = 0 at s and s
u = v = w = M = 0 at s I ands2
This results in eight constraint equations which when ex-
pressed in terms of power series solutions for eight sub-
domains along the meridian become
10 1
I a.s = 0 at s = sI and s = s2
i=o
73
10 -
Z b.s = 0 at s = s and s = s2
i=o
12
Sci.s = 0 at s = s1 and s = s 2
1=o
12 2 2 i2
E c.[(i)(i-1) + v(-n /sin 2 + i)]si- = 0
1=o
at s = s
and s = s2.
APPENDIX II
APPENDIX II
CONVERGENCE
In the use of the subdomain method, convergence
in the mean is desired. The main influence on convergence
is the choice of trial functions. For assured convergence,
the trial functions must be complete and linearly indepen-
dent [49]. The completeness property of a set of functions
insures that we can represent the exact solution provided
enough terms are used. The power series trial functions
used in this study are complete, for example, so that any
continuous function can be expanded in terms of them.
To demonstrate the rate of convergence for the
problem under study, a series of computer runs were con-
ducted to determine the dynamic perturbation response for
different numbers of subdomains. Dynamic stability response
for eight, ten, and sixteen subdomains was investigated.
The results are shown in Figure 11. The ordinate designa-
tion implies percent differences with respect to the best
approximation.
76
8
I
I
6
.- 4
4P
uJ
C
S-,
2
0 \
5 10
Number of Subdomains
Figure 11. Convergence Properties of Solution
15
APPENDIX III
APPENDIX III
FLOW DIAGRAM OF THE COMPUTER PROGRAM
For a better understanding of the computer program,
a key which may be helpful in going from the theory to the
program is given below:
FORTRAN Name Theory Description
N m Number of subdomains
NW n Circumferential wave number
Q Q Total number of n
AP a Semivertex angle of cone
RO p Material density
E E Young's modulus
V V Poisson's ratio
F s2 Conical slant height
BTT B Boundary condition matrix
STT [k] Linear stiffness matrix
DND [m] Mass matrix
PV {p} Load vector
PLUG {NL} Nonlinear stiffness vectors
X T Nondimensional time
FORTRAN Name
P
PO
T2
r
Nondimensional displacements
Axial load
Pressure level of Heaviside pulse
Upper time limit
Error bound
Online core storage matrices
Hamming's predictor-corrector
integration subroutine
Gauss-elimination subroutine
A flow diagram of the computer program is illustrated
in Figure 12.
U,V,Y
P
PSD
PRMT(2)
PRMT(4)
ORAN(I)
DHPCG
DGELG
Theory
Descriptions
80
VC
42'
V)
. o+
-- O--
(-) C 0
PS.-
LL0
L - - 0 \.
0
(-
- Y^ A4-
0
4- 00
0
=-3 4-
*U-
*I
u-
APPENDIX IV
APPENDIX IV
FORTRAN IV SOURCE PROGRAM
The computer program developed to provide the
numerical solution to this dissertation problem is listed
below. The program requires 182,000 bytes computer core
storage space. A typical nonlinear dynamic stability run
executes in approximately in seven minutes. The compile
time is three seconds.
00
00
00
00
00
Q00000
o0000
00000
S 00 0
0 0 *
N eo
. W w
%. Ww, .
0 f 0
0 I 0 I
0 m *.
-Q CO
. W. c.a 0
-SL 000
o *0
S-
N **Q
c C0 P,4
- 0W0 0
<
CL -
OM cU %0
o .oo
0> 0 00
9090
0 00
0000
0000
0000
0000
0000
cn -r n o
90QOO0000000909O0O
00 0000009 09009000
000000000000000000
000000000000000000
N N
p.- a
-*
c- <
0 N
~ ~ m N
Nf W
0 a
r~~~i~ N QC N0
-4 N
0 -
0.4 0.G 0b Y
a -4 Nm
W W* -M (\ '
00 0 00 ^ il a4
M ( * 0* N
N M *.
N N a a 0m N P*
Z)(41- ;;N 0.
- C). m NI :
X l a. ON w
V) *o *( Os ) 00 0- -00 G-
N 0 N N. 0.aaOO M.Jm 00 a
- - M VON-4-4 M V-4 (N
( N m Qo'-i4 * a NNE
0 0 0 *r- *4 *'Q4. ** 0
0. 0. 'T* 0.W W U) Nir o (\i\'4
( S - M U I M 0- 01
Nr N 4 0 -N NN 0 "- r ANNO W
0. NEMEm1mU N NN'
N N 4 N---> U
-0 00 0- 000@0-4) 0N OO -
- N N- N. W ** a- b- *0. *0 -
0' UON 0o 0 '--NO O m woo N ).
- *V) *tI) *v) 'o-M4, -- E
x-~ -- a'~^*a* W- o **M-. o;
WWWWW WW~e *WWW W* NWWWWL
* *a a** ** -W -MN N --4 9a *
x XNNXNNN4.4-Nr-. P-4 x1 NN- oX-
bO.-4 .4 b." .4 -.. 4 N -4 N N -4 N -N N- -
~---4 N~.r()m~ fnm(_)r.-P-. cu) >- a9~ ) o(1-
.. ..IAV) "***1 00 .>
~ZZZZZZ Z ON0ZZ 24ZZZ
000000000000^-0 OO 0OOOOQO
9-. *4 4 S... a a-- S-
00000 00 o0LUfJOCO.OQ mo O 000U0
N .-4 <4
-
00000
0090
00000
00000
en 0) (
00ooo00090
00000000
0o0ooo0o
Q00000000
900090009
990000000
000000000
Q00O000009
0oo00000o000
090000000
p 04 NMm -4rnQ %- W
in ttn o6 Lnin 6
oQoooooooo
oooWQooooo
ooooooooo
0000
0000
0000
0000
000000
000000
000000
m 0 0 r Q
10 a%0%QQ
0.
0%
m~W
OD 0
(A 0
N 00 0. o
N Z
0 N 0 ix
9 ** .S *
Lop) 00 00 o 0
QN Qe ** * *
*0~~ Z? OD m >bm
*0 t*< U & E 0 *00 0
NU > M 0 0 3c 3 X O
(M (1 . 0 Q** '9 -
C o- M P %- 0 It
&) Z
MN %Twm ob 0. X *
0- Lt. 'DNoc
1.- 0 9. 0
N 9. 3 x CO 0(o.o. -
0 * NNN p 1 e a .K .
- mN '0 9. 9. <\ Q 9 o
0ON CNo. ft f 0. z
V) < Cir .4 L m <_? >* ;nm <[
0. **QN U 0 Nc 0X 0M oc
oN U <- -s ) 0. D A 39 A. 0
M N Mtn 0 D OD 3
co N U L *>00 0 x.<
NO o N i r- 0 wI- I-N 0
** 00 m 000fn is it M 0 00 00 -
09 v. PQ.) (A) rQm m 0 -L) Do
a) N N-f 1- ft ;D u. u a% L- x. N ao N
en N (a a -4 -AN C.. U.. "0. cr.( u op'j. ol-9 w ) z
I I.- I.- ON u % w I Q. LU .C 0o <
V) V) V) # Q s-S inLLUfMV Q D -- CO ODPs CL
a* *N. s. 0.. U00' U0 V. Ua. -j -A 0 U- N
V) M mI- w 9. .4(A c-i Q.LI C- L u G- 0-4 *9. 6, CC) xo OD f
o. 4 --- 1---- 41 Q *Q '4)0 L X L.U O 4 9. *S u z x -
9.0 ,d /7 A- FN4 4.' u o) 0 .-dUeL (A 9( Z.~ VI D N0C 9O 9V
V) 0--//?\ LO V)V -L U m .u QC .
)( 9. NM 9. 9PSJ (ANs LA..LL0 9 (-3 )- 9. (M.a '0 -'-' 0
z .4 z z. J O C.C z z z z z e-) zzzzz -4- Aow a 9 w .
U W..4i0a W- 9.0 9u.0 9.LLO *uxI- U&'-4-U u 9 x%-.
m ( 4 *-4 b-4O b-d o-4 "4 s- U w o- o.i -4 "A "- "0 1-' 0% o w r- V
m) N -V AV )V nV )V iV V )t nV )V )V )V ) UJ-. ,- UJ m
(A UU O UluU U UOJU UWU^--4--WW*
0 WM 1- i-i -4 W 0-4-W- -P --i-W4- b. 4 OW i 6-4 -b4 U rl- L P- 00U
(N r> UJ < i U< W W tu L ) UJ UJ u ( t/ t ^ U< l < ui Uj LU UJ uL
1-1 u 6 .j -i b-i i-i j-i j- j b j-i i-ii -i >-< i-i i-i >i > -i i- -ii i^ i i-i j j i-j >i W~ (_ o W ou r t u
Uv- *> Q<)UUUUOUUUOUUUU(-UUUUUUU
(M UJW 0 0 00 0 0 0 0 000 0 0 0 00
00 00000000000000000000LJ0000(.30300
U N-4 -E
00
00
00
09
oo
00
OOOOQOOOo
Q000000000
000000900
r- f- ouor- o
00090 00900090000000
9 9009 090090090 00900
0000 000000000000000
~O~0'Q~0'0'0'Q~099
oQQ0000j~'~
0000000000000000000
909009090000 0000000
4
x
0
0
*
0 *
4 O
,0 0
0 0
7 0
< *o
-0o -
Z-. z
4m <
P0 (%
04
C% <(
(* *
< *
&-
*o -'
44n
S0
-. Q
0
(A
* 0
0
-- -
Q Q
000
- --L .
ac > Ic
o4 o
0. 0 o
1-
0L
Z 00
0 Lo QC
NrI O N c
.M I Q 00D
We .0 0 P-4 -t
uu x NM I" "- co
D-400 0 6 0
EoI OcO Q* *
f-, *,,, O e. e -
> t# I 'ft IM i of M -.0 -
4< U II Cn IIn II X X
MM M A -1,|, (N N N(
o09
00
oo
0 <
00
0 0
000
000
000
000
000
000
0 a
N 0
0 N
0 -)
1* 4
0 Qn c-
*~ -- 1-4 -
usOC W
M M1a. mm t- t-u -.l
N ** 1_^0U J
' *''0000La. X L.0
*O OM ) L
V W' ------- mo -u -
gq *S P: Li Ic u4 0 xx-4o
a:W(000000J004100
W L L. ***O. *U* U.n
-iMuAMMM>** -MNXXXXXX
O(x QacxxQQQooooo
N\ 0 -1-4 Nr %C '~ w \.'~ <^
P- r6-t- P-I-w- -I- - -
44 M qX
fn
OQOOOQQOOQO
o OQQ.OOQ~QOD
V0 44040.0 W 0Q 0-4 -4 -4
9O00000O0000000000000000
O000000 OO0000000000 000
0000000000000000oo00000
0000000000000000000000
00000000000000000000003
Q 0 ~r0i ca Q D.JQOQQ
m 04 .-'VK40 WMMNr4 9-4 .00
N ON --0 ONa' 0440C
N ~ ~ ~ w ~ o0 0% N0 r0 %N
0n o r- u- 0q Q Goi 0 0 n M Q Q N00 0 00 + (
000000009000000000+--
~~~~C gg 0 40 ai 0 0 Q 0oo 0 Q Q N -00 W00 0
-- - Q'.4 P4 0 .4 C- 0 0 0 )C ;l + 4 N
w i 1 f% a 0 ~N fn W Ln W *, I 0 1 *41 M w 4V) Z V)zeo
a~-~ ~ It u It of t 0 1 11 It of 1 to HQ
11z of
w 1 -4
x a,
0000000000
0000000000
00000000000
0000000000
Q m -d N M w in Q r,
w w T oow w ro-
V4 4 ooo 04 4o .4Po
0000
0000
0000
0 to
0000
00000000000
00000090000
oooooooooo
0000000000
Nomoronooo CPoroC1oo4
00009000
00090999
00000000
0.4 W4 04 0104 -.4 0" 0.4
00000000
00000000
"4-0
V PV)
VPZZ
1/) Z Z
ZNO
Z N 0
..v a (
9M ..I 9.
N**
*V)
-M a- -<
0 ~~ fM* l/
0 OLLJO
0 (^ 0 f
G >
-< t ri 2L
I r-4 ;
, I- -I --
'-0 ftN (M
Lrr'
*q (f) L40U1)
0 L400
*0 N(/)V
4 0 (MO7tJO
< U0 U
dO- w--
II 1t > > II
- If )0.
-=WNO.ZQ
.- ) A. x
* ftZNVi/a.VP-r
C-4M ev || If
-- 04 _1J _j C
S-4 Z /) .J
0 t0j di
V0
ZV
V*
V)
z
0
z
0z
CM4
tC)
9.
0
o.
N
0C
*-
it I II
L)J V) V) vJ LO
< V) V) < 00(f
0 .4 4 a
V)
,)
9.1
0z
. 0
I
0.
z
* -
SM
- 0
0 0
f *9
- -
I 0
'-4 *t 0
VP U ff
. M ,- -
-J NN
..) t
Z O Z Z
S. I-
aM
0 N ft
0 ) 0.
(a UN M
I Q- *
-0 6NC
N Q0 Q
It of >- > >.
--) N It M
D Ow D X X
u- 0 N N -l < M II
t -MCu uv z Q uJ-lQ
0
0
0
"4
0-
*
M-
I t
0-
0 *
0p4
"4
t
i- n
l -
1-1
0 0
0
0 .N
II -
*- Z z.
ut 0 of
RSS
