Recent advances in reduction methods for instability analysis of structures

Comp~ter5 & Str~cture., Vol 16. No. I--4, pp. 67-81L 1993 Pnnted in Great Britain

004%79491831010067-14503.00]0 Pergamon Press Ltd.

RECENT ADVANCES IN REDUCTION METHODS FOR INSTABILITY A N A L Y S I S OF STRUCTURES AHMED K. NOOR# and JEANNE M. PETERS~ George Washington University Center at NASA-Langley Research Center, Hampton, V A 23665, U.S.A. Abstract--Status and some recent advances in the application of reduction methods to instability analysis of structures are summarized. The aspects of reduction methods discussed herein include: (a) multiple-parameter reduction methods for instability analysis of structures subjected to combined loads: (b) use of mixed models in conjunction with reduction methods for instability analysis: and (c) application of reduction methods to instability problems with displacement-dependent and nonconservative loadings. Numerical examples are presented to demonstrate the effectiveness of reduction methods in instability analysis. Also, research areas which have high potential for application of reduction methods are identified. radius of curvature number of basis vectors (reduced degrees of freedom) current stiffness parameter defined in eqn (18) [s] generalized stiffness matrix defined in eqns 06) U total strain energy of the structure Ul. U2, W displacement components in the coordinate directions {xt vector of nodal dispiacements )fl, X2, "~3 Cartesian coordinate system Ix3 normal to the middle plane of the plate)

NOTATION

A cross sectional area d typical design variable E,, Er elastic moduli of the individual layers of a laminated plate in the direction of the fibers and normal to it, respectively [F] linear global flexibility matrix defined in eqns (16) / rise of shell (or geodesic dome)

(/(X, a,, z,)},

{f($, Ai, A:)} vectors defined in eqns (1) and (6)

(/.(/-/. X, a,, ag}, {/x(H, X, x~, a.,)} vectors defined in eqtas (16) GLr, Grr shear moduli in the plane of the fibers and

{ Gfx)}, { ~(x)}

h I

[K1 [K], [~ [K"q, [K ~:)]

[r], [~

[rx], [r.], matrices of basis vectors defined in eqns

normal to it, respectively vectors of nonlinear terms defined in eqns (1) and (16) vectors of nonlinear terms of the reduced systems defined in eqns (6) and (21) thickness of plate, shell or ring moment of inertia of the cross section global linear stiffness matrix of the structure linear stiffness matrices of the reduced systems defined in eqns (6) and (21) geometric stiffness matrices defined in eqns

(5), 114) and (20) :t. hi, A: path parameters VLT major Poisson's ratio of the individual layers of a laminated plate rotation components of the middle plane of the plate vectors of undenermined coefficients of the reduced equations, defined in eqns (6) and (21)

(2) [K'"q normalized load stiffness matrix defined in eqns (26) [R'"q load stiffness matrix of the reduced system L

{~(H. X)} N~, N,.~

NI¢, n

{Q},{¢"},{¢:,} { O}. { Qm}, { O,.~} q, q,, q~

I. INTRODUCTION Since the development of the general theory of elastic buckling and postbuckling by Koiter in 194511], considerable progress has been made in modeling various instability phenomena in structures. Review of the many contributions on this subject is given in a number of survey papers (see, e.g. [2-8]). Although several finite element programs are now available for predicting the critical loads of practical structures (see, e.g. [9]), the establishment of reliable and efficient finite element techniques for instability analysis of structures with complex geometry and/or loading remains a challenging task and continues to be the subject of intense research effort at the present time. Most of the finite element applications to instability analysis of structures are based on using either (a) incremental/iterative approach, or (b) asymptotic (or perturbation) method (see, e.g. [10, 11]). The latter method is an adaptation of Koiter's perturbation procedure for the study of the immediate postbuckling response and the sensitivity of the structure to fabricational imperfections. The algorithmic tools for finite element instability analysis encompass three distinct aspects: ta) nonlinear (or

defined in eqns (30) side length of the plate or planform dimension of Geodesic dome vector of nonlinear terms defined in eqns (16) in-plane normal and shear stress resultants acting on the plate critical value of N~ for the case of pure compressive loading (NI: = 0) total number of displacement degrees of freedom in the discretized structure normalized external load vectors normalized load vectors of the reduced system load parameters associated with {Q}, {Qm},

{O':'}

critical value of q critical combination of q~ and q~ tProfessor of Engineering and Applied Science. ~Programmer Analyst. 67

~8

', K. Noo~ ami J. M. Pr~rF~ ~

linearJ analysis of the prebuckling state. (b) determination of the critical points (e.g. bifurcation and limit points) on the equilibrium path, and (c) tracing the postcritical response. Special numerical algorithms have been proposed for overcoming the difficulties associated with commonly-used iterative techniques near critical points (see, e.g. [12-17]). In spite of all these developments, finite element instability analyses of structures with complex geometry and/or loading require excessive amounts of computer time even on present-day large computers and thus are very expensive. Considerable work is currently being devoted to reducing the computational effort required for instability analysis of large complex structures. This work includes: (a) using simple mathematical models that capture the major effects in the response of the structure; (b) using efficient numerical algorithms for solution of nonlinear equations and extraction of eigenvalues [1820]; (c) reducing the size of the analysis model by exploiting all the obvious and nonobvious symmetries exhibited by the response (see, e.g. [21]); and (d) reducing the number of degrees of freedom of the initial discretization through the use of reduction methods. The essence of reduction methods is to limit the deformation modes of the discretized structure to few known modes (basis vectors or global Rayleigh-Ritz approximation functions). Due to the high potential of reduction methods for nonlinear analysis increasing interest has been shown in the application of these techniques to nonlinear static and dynamic problems[22-30]. This paper reports on some recent developments of reduction methods for instability analysis. Discussion focuses on a number of aspects of reduction methods including: (a) multiple-parameter reduction methods for instability analysis of structures subjected to combined loads; (b) use of mixed models in conjunction with reduction methods for instability analysis; and (c) application of reduction methods to instability problems with displacement-dependent and nonconservative loadings. Numerical examples are presented to demonstrate the effectiveness of reduction methods in instability analysis. Also, research areas which have high potential for application of reduction methods are identified. 2. MULTIPLE-PARAMETER REDUCTION METHOD8 FOR STRUCTURES SUBJECTED TO COMBINED LOADS

In many practical problems, the structure may be subjected to various combinations of loads, and multiple independent loading parameters may be required to describe the external loading. In this section the application of multiple-parameter reduction methods to instability analysis of structures subjected to quasi-static loading is discussed. For simplicity the external loading is assumed to be conservative and is represented by two independent parameters q, and q2 only. The cases of displacement-dependent and nonconservative loadings are discussed in a subsequent section. A total Lagranglan formulation is used and the structure is discretized by using displacement finite element models. The advantages of using mixed finite element models, with the fundamental unknowns consisting of both stress and displacement parameters are outlined in the succeeding section.

2.1 Governing finite element equations The governing finite element equations of the discretized structure can be cast in the following form:

if(X.,~,,,~_.li=[K]{X}-GiX}-- ;,iQ' - ~ <' :})

where [K] is the n by n linear global stiffness matrix; a is the total number of displacement degrees of freedom: iX} is the vector of nodal displacements; {G(X)} is the vector of nonlinear terms: Q'~'~ and {Q'->, are nor~ realized load vectors; q, and q~ are independent load parameters; and ~., and .<, are independent path parameters which may be identified with loading, displacement or arc-length parameters in the solution space. In the latter case, eqns (1) have to be augmented by two constraint equations relating ~a,,a,) to both ~.V~ and (ql, q2); see [29]. The deformation of the discretized structure is described by the set of n displacement variables, components of the vector {X}. An equilibrium state corresponding to given values of q, and q: is associated with a point in the in + 2) dimensional load-displacement space spanned by the components of the vector IX} and the load parameters q, and q2. The totality of the equilibrium points constitutes a two-dimensional equilibrium surface. In contrast, if the loading is described by a single parameter, then the totality of the equilibrium points constitutes a curve in the solution space. The standard approach for the solution of eqns i i) is to fix the value of one of the two parameters q, or q_, and vary the other: or to choose a functional relationship between q~ and q2 which is dependent on a single parameter q. In either case, the entire solution corresponding to the chosen combination of q, and qz (which is effectively dependent on a single parameter~ constitutes a curve on the equilibrium surface of the structure. 2.2 Identification and determination of stability boun-

dary For certain combinations of the load parameters q,, q: an instability occurs and the corresponding point on the equilibrium surface is called a critical point. The totality of the critical points in the load space constitutes the stability boundary, which separates regions of stability and instability. The stability boundary defines an implicit relationship (or an interaction diagram) between the two load parameters q, and q> For conservative loading which can be described by a single load parameter the critical points on the equilibrium surface are either limit or bifurcation points, provided the structure is elastic and has no imperfections. For imperfect structures, bifurcation points degenerate into limit points. Bifurcation and limit points cannot describe adequately the instability of the structure in the presence of more than one independent load parameter, and new classification of critical points into special and general has been suggested in [31]. The former is a genuine bifurcation point. The latter is a limit point but under some conditions it appears as a point of bifurcation. Critical points are detected by the vanishing of the determinant of the tangent stiffness matrix [[K] + (OGIIOX~)] where [ and J have the range 1 to n. In the determination of the critical points two cases are to be distinguished: (a) When the deformation in the prescribed state is either rotation free or the effect of rotations can be neglected in the prebuckling analysis. Examples of this situation are provided by flat plates subjected to in-plane loading and circular rings subjected to pressure loading; (b) When the rotations in the pre-

Recent advances in reduction methods for instability analysi~ of structures buckling state are important, such as in a spherical cap subjected to concentrated load. Case I. Prebuckling rotations are unimportant, The stability boundary can be obtained by linearization of the governing finite element equations, eqns (1), about the initial undeformed configuration ql = q2 = O, {X} = 0 and solving a linear algebraic eigenvalue problem of the form:

[IK] + 0,i "'l +

0

where [F] is an n by r transformation matrix whose columns represent global functions or basis vectors and {0},., is a vector of unknown coefficients (reduced degrees of freedom) which represent amplitudes of displacement modes. The matrix [F] is given by:

'

a2X laAl'J

_

""

' (5)

(2)

where (c~, 4,_) represents a critical combination of the load parameters q,, qz and {A'} is the associated eigenmode. The geometric stiffness matrices [K"] and [K("q are obtained by setting {X} equal to {aXlaq~}o and {OXIOq,.}o. respectively in [aG~/aXj, where the subscript 0 refers to the value of the vectors at q~ = q., = 0. The standard computational approach for determination of the stability boundary involves: ta) Evaluation of {OXlaq~}o and {aXlaqJo through solution of the linear equations:

f ( aX ]

6Q

The equations used in evaluating the basis vectors are obtained by successive differentiation of the governing finite element equations of the initial discretization, eqns (1). For the single-parameter reduction method the explicit form of these equations is given in [25] and is not reproduced herein. The Rayleigh-Ritz technique is then used to approximate eqns (1) by the following reduced system of nonlinear algebraic equations in the unknowns {tO}: ~(6, a,, ~,)} = {R]{6} • {d(q,)} - q,{Q'"} - q~{O'"'}

(6) aX where

(b) Generation of the geometric stiffness matrices [K"'] and [K(:}]. (c) Repeated solution of the linear generalized eigenvalue problem given by eqns (2) for different combinations of (~t and (~.~.Obviously, if an accurate determination of the stability boundary is desired, eqns (2) will be solved many t~mes at a considerable cost, even if an efficient eigenvalue extraction technique is used. Case 2. Prebuckling rotations are important. The stability boundary can be obtained by using either one of the following two approaches[32]: (a) Gradual incrementation of the load parameters (q~, q,_). solution of the corresponding nonlinear equations, eqns (1), and evaluation of the determinant of the tangent stiffness matrix [[K] + (OGdOXs)]. The critical points correspond to a zero value of the determinant, or (b) Finding an equilibrium configuration near the stability boundary, and assuming the incremental displacements from that configuration to be proportional to the load increments. This leads to a quadratic eigenvalue problem in q~, q:. Efficient numerical procedures exist for solution of quadratic eigenvalue problems. However, the spurious complex roots which may be present near the real root can result in slowing the convergence of iterative techniques. If the reference state is sufficiently close to the stability boundary, a linearized eigenvalue problem, similar to eqns (2), can be used for determining the stability boundary. The vectors {aX/aq~} and {aXlaqJ are evaluated at the reference state and (4~, ~:) are measured from that state (see, e.g. [33]). 2.3 Basis reduction and reduced system of equations The nonlinear solution {X} in eqns (I) is approximated over a range of values of the two load parameters ql and q:, by a linear combination of a nonlinear solution {X} for a particular combination (ql, q j and a number of the path derivatives of {X} evaluated for the same combination (q,, qJ. The approximation can be expressed by the following transformation: {X},,., = [FL.,{~},.,:

r~ n

(4)

[/(] = [F]r[K][F]

(7)

{d(qn} = [rf- {G(O)}

(8)

{O(''} = [F]T{Q'''}

(9)

and {G(@)} is obtained from {G(X)} by replacing {X} by its expression in terms of (@}. If the prebuckling rotations are unimportant, the stability boundary is determined by first generating initial basis vectors at q~ = q.~= 0, {X} = 0 and then using Rayleigh-Ritz technique to approximate the linear eigenvalue problem, eqns (2), by the following reduced system of equations:

where [ / ( " l = [~]T [K"'] [~

(12)

[£,2,] = [~T[K,: q[p]

(13)

and [r'] is the matrix of initial basis vectors defined as follows:

t aa,-Jo t aa---~do aA, /e,2--~-Jo"'' ] (14) and subscript 0 refers to the value at q~ = q: = O. Once the reduced eigenvector is obtained, the buckling mode for the full system is generated by using the following equations:

{.~'}= [h{6}.

(Is)

3. USE OF REDUCTION METBODS IN CONJUNCTION WITH MIXED MODELS

In a number of applications of reduction methods with

70

~,. K. NnnR and J. M. PETERS

displacement finite element models it was found that the accuracy of the stresses predicted by the approximate reduced equations was, in general, lower than that of the displacements. An effective way to remedy this drawback is to use reduction methods in conjunction with mixed (multifield) finite dement models with the fundamental unknowns in the initial discretization consisting of both stresses and displacements. The basic idea of this approach as applied to instability analysis of structures subjected to combined loading with multiple independent load parameters is discussed subsequently.

discussed in the preceding sections, the basis vectors are chosen to consist of a nonlinear solution and its variousorder derivatives with respect to the path parameters a,, ,L as follows: {H}= [~;J/O}

~19,

where 0

H

a

H

3.1 Governing finite element equations A total Lagrangian formulation is used and the structure is discretized by using a two-field mixed model with the fundamental unknowns consisting of both stress and displacement parameters. For simplicity, the external loading is assumed to be conservative and is represented by the two independent parameters q, and q> The governing finite element equations consist of both the constitutive relations and equilibrium equations and can be cast in the following form:

[ [S] F

[&(H, X, a,,

(20) The recursion formulas used in evaluating the basis vectors are obtained by successive differentiation of the finite element equations of the initial discretization, eqns (16). A Rayleigh-Ritz technique is used to replace the governing finite element equations, eqns (16), by the following reduced system of equations:

II.~(H, X)jl

{/(~, x,, a~) = [~q {~} + { q~(O)}-

q,{O.'"}- q~{O.'~>}

(21)

where where [F] and [S] are the linear global flexibility and generalized stiffness matrices; {H} and {X} are the vectors of stress parameters and nodal displacements; {~(X)} and {A/(H,X)} are vectors of nonlinear contributions. The other terms are defined in connection with eqns (1). Note that if the control parameters .~ and X2 are selected to be generalized arc lengths in the solution space, then eqns (16) have to be augmented by two constraint equations relating (,L, a,) to {H}, {X}, qt and q2 (see [29]). If the prebuckling rotations are unimportant, the stability boundary can be obtained by solving the following linear algebraic eigenvalue problem:

i,=ljt{ }=o

iS] T

(17) where [K <''] and [K~>] are the geometric stiffness matrices obtained by setting {H}, {X} equal to {dHIOq,}o, {a,,~0qt}o and {aXldqz}o, respectively in [o./ddOXj]; and {H/X} is the ei~enmode associated with the critical combination (q!, q2). If, on the other hand, the prebuckling rotations are important, the procedure outlined in the preceding section is used and the critical points are obtained by finding the points at which the current stiffness parameters[34] is either zero or has a discontinuity, where

1

OX

(18)

a = 1 or 2 depending on whether q~ or q_, is incremented, and So is the value of S at qo = 0. 3.2 Basis vectors and reduced system o[ equations As for the case of displacement finite element models

[~,-] = _ [FHlr[F][Fu] + [FHlr[S][F×] + [rxIr[S]r[r~] (22)

{<,3(~)}= [F.IT{~4)} + [rxl r(~a(4)}

(23)

{(~m} = [Fxlr{O,,}

(24)

{0'-"} = [rx]TfO'~}.

(2~)

and

In eqns (23), {W(~)} and {.~(~b)}are obtained from {~X)} and {dt(H,X)} by replacing {H} and {X} by their expressions in terms of {$}, eqns (19). The reduced eigenvalue problem corresponding to eqns (17) has the same form as eqns (11). 3.3 Assets and liabilities of using reduction methods in

con]unction with mixed models If the proposed approach for applying the reduction method in conjunction with mixed models is contrasted with the corresponding displacement approach, the following major advantages can be identified. (1) Simplicity of generation of basis vectors. For geometrically nonlinear problems, the nonlinear terms in the finite element equations of the mixed method have simple mathematical structures, and are bilinear (or quadratic) in the nodal parameters[35,36]. In contrast, the nonlinear terms in the displacement approach are cubic in the nodal displacement parameters. As a consequence of this, the evaluation of the path derivatives for mixed finite element models involves less arithmetic operations than that of the corresponding displacement models. (2) Better approximation properties. Numerical results reported in Refs. [28,33,35] for the single-parameter reduction method have shown that for a given number of basis vectors, the accuracy of the solutions obtained by

71

Recent advances in reduction methods for instability'analysis of structures the reduction method-mixed model is higher than that of the corresponding displacement approach. Thxs can be attributed to the explicit approximation of {H} in the mixed equations (through the use of [F,]). While there are a number of advantages in using the reduction method in conjunction with mixed models, two major difficulties also arise: (a) the first results from the large number of degrees of freedom used in the mixed method. This leads to a substantial increase in the number of simultaneous equations used in generating and updating the basis vectors, and (b) the second difficulty is due to the nondefiniteness of the matrix of the algebraic equations used in generating the basis vectors. The first difficulty can be alleviated by using mixed models with discontinuous stress fields at interelement boundaries. Such mixed models were shown in Ref. [36] to be equivalent to reduced/selective integration displacement models and have better performance than mixed models with continuous fields. Moreover, the discontinuous stress field and its path derivatives can be eliminated on the element level, thereby considerably reducing the size of the system of equations used for evaluating the path derivatives. Also, the use of mixed models with only one continuous field at interelement boundaries simplifies the implementation of the reduction method in existing nonlinear programs based on the displacement formulation. Note that for pin-jointed truss structures the global flexibility matrix [F] in eqns (16) is diagonal, and therefore, the elimination of the forces and their path derivatives, in the process of evaluating the path derivatives of the nodal displacements, can be easily done. 4. DISPLACEMENT-DEPENDENT AND NONCONSERVATIVE LOADINGS

In recent years finite element instability analysis of structures subjected to displacement-dependent loads has attracted a great deal of attention. Among the many publications on the subject, mention may be made of Refs. [37-41]. In all the cited references displacementdependent loads are each defined by an a priori assumption that the load acting on the structure is calculable from an assigned function of at most the generalized displacement of that point and its neighborhood. The size of the displacement is unrestricted. However, time and time derivatives are excluded from such a function. Time-dependent loads are by nature nonconservative, and a stability analysis for structures subjected to these loads can only be made on the basis of equations of motion[42]. A general treatment of the mathematical properties of pseudo-static displacement-dependent !oadings is given in [43], wherein criteria for conservative loading are derived and applied. Consideration of displacement-dependence of the loading may result in significant changes in the critical values of the loads for certain types of structures (e.g. circular ring subjected to hydrostatic pressure). Within the framework of the finite element method, the displacement-dependence of the loading (or follower-load effect) is reflected in the so-called load stiffness matrix. This is a preload stiffness. It is nonzero only if the force is nonzero and is, therefore, relevant to bifurcation buckling and nonlinear analysis. The load stiffness matrix is symmetric for conservative loads and is unsymmetric for nonconservative loads. An example of the latter case is provided by shells with loaded free edges[41]. The load stiffness matrices associated with pressure ioading,

centrifugal load. frictional drag and cable toad are given in [38]. The admissibility of symmetrizing the unsymmetric load stiffness matrix in the instability analysis of structures subjected to nonconservative loadings has been studied in [40,411. In this section the application of reduction methods to instability problems of structures subjected to pseudostatic displacement-dependent and nonconservative loading is discussed. For convenience, the external loading is assumed to be represented by a single load parameter q; and as in the preceding sections a total Lagrangian formulation is used and the structure is discretized by using displacement finite element models. 4.1 Governing finite element equations The governing finite element equations can be cast in the following form: F

"1

k

J

{/(x. ~)t = / [ ~ : 1 - q[K'")]/{X} + {G(X)}- q {Q} = o (26)

where [K ~"'] is a normalized load stiffness matrix (independent of the load parameter q): {Q} is the normalized load vector; and )t is a path parameter. The other terms are defined in connection with eqns (I). Note that for the case of nonconservative loading the matrix [K'"q is unsymmetric. If the prebuckling rotations are negligible, the bifurcation buckling load can be obtained by solving the linear eigenvalue problem defined by the equations: [[K]+ 0([K'"] - [K'q'])] I.~} = 0

(27)

where q is the critical value of the load parameter q; {X'} is the associated eigenmode; and [K'"] is the geometric stiffness matrix obtained by setting {X} equal to {OX/Oq},, in [oGdOXj]. 4.2 Reduced system of equations As in the case of multiple-parameter reduction method, the basis reduction is accomplished by using the following transformation: {J~'}= [F]{qA.

(28)

If the basis vectors are generated at q =,~ =0. the matrix [I'] is given by:

fo=Xl___,_Io~Xl~

[fox/__ jr] = 1.l o a / o / o ~ - / , , / o A - / .

]

"

(29)

The recursion formulas used in evaluating the basis vectors are obtained by successive differentiation of the finite element equations of the initial discretization, eqns (26). Since for q = 0, the matrix [K 'q'] does not appear on the left hand side, the generation of all the basis vectors involves the decomposition of the symmetric global stiffness matrix [K], The reduced equations corresponding to eqns (26) are given by:

= [telL

J

q 0/= 0 (30)

:"~

;',

K. NOORLindJ. 'vI. PETERS

where

P

[R'>] = [r] r [K"][F].

3 I)

The other terms are defined in connection with eqns (6). For nonconservative loading the matrix [/~q)] and hence the reduced equations, eqns (30), are unsymmetric. However, due to the small number of these equations, no approximation or symmetrization need to be made in the process of their solution. The reduced eigenvalue problem corresponding to eqns (27), has the following form:

E = 2.%5 × ].09 ~'4,;m2

=0.33

v

R = 0.254 ,n. a = 6.35 x 10 ;2 m, h

=

6.35 x 10 4 m.

c( = 14.4780 0

Once the reduced eigenvector is obtained, the buckling mode for the full system is generated using eqns (28).

Fig. 1. Shallow spherical shell used in the present study.

5. NUMERICALSTI0"DIF.,S Numerical studies which demonstrate the effectiveness of the reduction methods presented in the preceding sections for the instability analysis of structures have been presented in [28-30]. Herein the results of four

typical problems are discussed. The four problems are: (a) instability of shallow spherical cap subjected to combined external pressure and concentrated loads; (b) bifurcation buckling and postbuckling of laminated

I

Full system(148D.O.F.) 0 +

14 vectors I Reduction 20 vectors method Ilil

I

i III

p

o o.,,.,o-3

l

q1=oA,,ir ~ '

I

"

_~00000

q2 = FR 1.5 Eh3 1.0

x 3

o.s

-"~,

~O 0

j J ~ ~ , ~ ~ ~ ~ m 0

-0.1

/~1-o t t_ |

/

Wc/f

-0.2

1831 x 10.3 ,.0.3

P

" _,.,,oOOOOOOOoo

q2 = ~

r- 611= 0.915x 10.3

Sh30.51"0If ( ~ ~ O C / oI~0 0,0I ~or~-~,/-ql=1"373Xl°'3 ~0~ O ~ ~z ~r q 0

2.5

5,0 UR Eh4

7.5

1 : 1.831x 10-3 10.0

Fig. 2. Accuracyof solutions obtained using multiple,parameter reduction method. Shallow spherical shell shown in Fig. 1, ¢ = (pRIEh),

Recent advances in reduction methods for instability analysis of structures composite plate subjected to combined in-plane compressive and shear loadings; (c) instability of a geodesic dome with two independent loading parameters; and (d) bifurcation buckling of a circular ring and a semicircular arch subjected to hydrostatic pressure loading. The first two problems were selected to test the effectiveness of the multiple-parameter reduction method when applied to structures subjected to combined loadings. The third problem is used to assess the merits of using reduction methods in conjunction with the mixed formulation, and the fourth problem gives an indication of the potential of reduction methods when applied to problems with displacement-dependent loadings. 5.1 ShaUow spherical cap subjected to combined pressure and concentrated loads The first problem considered is that of the axisymmetric deformation of the clamped shallow spherical cap shown in Fig. I. The shell is subjected to combined pressure and concentrated load at the apex. Experimental and analytical results for the critical combinations of the two loads are given in [44]. 180[--

"~

Full system (14~ D O.F.i

I C

ga

9 vectc~rs t

R~it~t on

2 Er3

o

o7

J

I

14

7.1

~

! x 10-3 ~2.~

. 0R %--~-ff

Fig. 3. Accuracy of stability boundary obtained using multipleparameter reduction method. Shallow spherical shell shown in Fig. 1.

Due to axial symmetry of the shell and the deformation, only one meridian was considered and was modeled by using ten shear-flexible shell elements with quintic interpolation functions for each of the two displacements and rotation components (a total of 14~ nonzero degrees of freedom). Twenty basis vectors were first generated for the unloaded shell [q, = (pRIEh) = O, q~ = (PR/Eh ~-)= O, {X} = 0], and were thus obtained by solving a linear system of finite element equations. The basis vectors included all the derivatives of {X} with respect to q, and q2 up to the fifth order, inclusive. The vectors were orthonormalized using the Oram-Schmidt technique in order to improve the conditioning of the resulting reduced equations. The same set of basis vectors was used for determining the stability boundary and tracing the nonlinear axisymmetric response of the shell associated with various combinations of q~ and q2. An indication of the accuracy of the predictions of the multiple-parameter reduction method for five different values of q~ is given in Fig. 2. As can be seen from Fig. 2, the transverse displacement w~ and the total strain energy U obtained using fourteen basis vectors accurately predict the limit points but, for smaller values of q,, they deviate from the full system solutions in the post-limit-point range. On the other hand, the solutions obtained using twenty basis vectors are in close agreement with the full system solutions in the entire range of loading considered. The stability boundary of the shell was obtained by determining the value of the limit-point concentrated load parameter qz corresponding to various values of the external pressure load parameter qt. When the pressure load parameter q, is less than 0.34 × 10-~, no limit type behavior is exhibited by the shell for any value of the concentrated apex load. The accuracy of the stability boundary obtained by the multiple-parameter reduction method using 9. 14 and 20 basis vectors is shown in Fig. 3. As can be seen from this figure, the 9 and 14 basis vectors do not predict the limit point concentrated loads

L = 0.254 rn h = 2.114 x !0 -3 rn

cF:i I/;I,

Properties of in(livi0uat lae~rs E L = 13.10x lO10Ntm 2 E!. = 1.303x 1010N/m 2 GL.[ = 0.64] x 1010 N/m 2

GTT = 0.261 x I010 Nlm2 VLT : 0.3~ Fiber orientation:

(±451021T-451902)s

I----#---

Bou nde r'y. co nditions Eage x I = L/2:

w

= ~1 = ~)2 = 0

Ix3 w

At x I : x 2 : 0: u I : u2= ~i = ¢2= 0

E(lges x 2 = :t:L/2: w = ~1 = ~2 = 0 . At x I = L/2. x 2 :

u2

0: u 2 = 0

/

xiP'

Fig. 4. Sixteen-plygraphite-epoxyplate used in the present study. CA5 I~:~14- r

73

74

\. K, N(;OR :rod J. M. PFTERS

Table I. Accuracy and convergence of minimum eigenvatues ~ obtained using multiple-parameter reduction method Sixteen-ply graphite-epoxy plate subjected to combined compressive and shear loadings (see Fig. 4/ N1 "-7"-Nlcr

Eigerrvalues

Number of Ehasls Vectors

(Sylne ~ric M~des)

!

9

- !4

1

Full System (792 D ' O ' F "

20

NIL~ a) Pure compressive loading ql = - E_h 3 1.0

~i)

-37.921

4,(2)

-89.634

! ~37.912

-37.912

-37.912

-84. 702

-84,701

-~,9~

24.913

~4~13

i -84,S48

!

1 b) Combined loading o. 75

~l~

0.50

~_~.,60

!

~2)

-27.797

i -27.761

-27.761

-27.761

~I)

36.555

36,441

36,439

36.429

q~2)

-39.495

! -39. ,~67

-39.466

-39,466

I

!

olz5 ~2)

-48.886

45.6,

1 -48,864

-48.86~

I

-48.86., l

"~

c) Pure Sheer Loading q 2 " -''.'3 I 1 ETn O.

~i)_

153.7

54.172

3.726

53.726

,7.069

-37 .069

%

q2 -(2)

-57. 094

I -57,0

t:-~

i

associated with small values of the external pressure loading. On the other hand, the twenty basis vectors accurately predict the entire stability boundary, 5.2 Bifurcation buckling andpostbuckling of a laminated composite plate subjected to combined loading The second problem considered is that of the bifurI

Full system

•

,7~ 0.o.~

i

+

I

........

2)| ~l~llm)l:N

tlvgtorsi '~"~

- |

i [ "~-

-"~-..~

'~1~

[~,~-

-" ,I * ~ - ~ ~ - c ~ "

m ~ , ~ --

-,.~\

~F~z t-~ 01i ~

~

h~ = -"-W %~-

/

_

_

_~j~~~ -0zs

~

~...@i

~,-

075

,~

..,r~.

~o

~'~,, Fig. 5. Accuracy of stability boundary obtained usmg multipleparameter reduction method. Sixteen-ply graphite-epoxy plate shown in Fig. 4.

ca,v,, ~,,,,.,ding and postbuckling analysis of a 16-ply graphite-epoxy square plate subjected to combined edge compressive and shear loadings. The material and geometric characteristics of the plate are given in Fig: 4. Both the bifurcation buckling modes and the postbuckling responses of the plate exhibit inversion symmerry characterized by the following relations: u,(x,, x,_) = - u,~(- x,, - xH w(x,,x,.)= w ( - x ,

-x,_)

a = 1.2.

(33)

~.(x.x..)=-~.(-x,,-xq

The above symmetry relations were used in c o n junction with the procedure outlined in [21] to reduce the size of the analysis model to half the plate. A grid of 3 x 6 shear-~xible elements was used with bicubic Lagrangian interpolation functions for each of the generalized displacement components (a total of 792 nonzero degrees of freedom). Typical results are presented in Table 1 and in Figs. 5-7 and are discussed subsequently. The basis vectors were first generated for the unloaded plate (q, = q , = O , { X } = O ) , and were thus obtained by solving a linear system of finite element equations° Then the reduced stiffness and geometric stiffness matrices were generated and used to determine the stability boundary of the plate. An indication of the accuracy of the stability boundary obtained by the multiple parameter reduction method is given in Fig. 5, and the buckling mode shapes associated with different cornbinations of c], and ~-. are shown in Fig. 6. The con-

Recent advances in reduction methods fer instability analysis of structures --~

b

"~

b

~-

- - ~ ~ 3 ~ ~" ~ .

%

op - - - - - - - ~ -

:~,'/,;zz; ,

/~?~((s,,,/

(1) q2 " 53.726

(2) CI2 • 45.673

~5

)111,!11

~.~q,//,,'i

U ,

N1/Nlcr N1 • 0

N 1 ~ 0.25 N l c r

i i

_

.~

(3} ~l 2 = 3 6 . 4 3 9

N1 = 0.50 N1¢r

N12 -7,,

b

L E:

•

';

1

"J~:"

ii ,:t~il :

" "L" " "

: .) i -a(,l) ~I2 = 24.913 N I = 0.75 NIe r

.)

II (5) q~ - 0 N I • NI¢ r

16) ~12 - -39.466 N 1 • 0.5 NlC r

Fig. 6. Normalized contour plots of w for various buckling modes (associated with different combinations of compressive and shear Ioadings). Sixteen-ply graphite.epoxy plate shown in Fig. 4, t~, = (N~,L:/Erh~), N~, = - 37.912(Erh3/L2). vergence of the lowest two eigenvalues ~'~ and ~:~ with the increase in the number of basis vectors is shown in Table 1 for different combinations of compressive and shear loadings, including the cases of pure compressive loading (q,_ = 0), and pure shear loading (q~ = 0). As can be seen from Table 1, the errors in the lowest eigenvalues obtained by using 14 basis vectors (first, second,

80

d"/+4"/

N1/N1

•r

,.X,

|

0.50 ~

~\~*'~,y"

~+.-.~.'~,~,~ q2

Nt2 ["1'.. . . . . . . .

40.

O

1.0

2.0

: '

a.o

Wc/h N1/Nlcr ~0.25 Bending energy~ ~/~0.50

Bo

:

oo

r

. q2 40.

2o

0

~W ~

p

\

--o.so

LTotaUstam energ~,~O.715

~

I

/P("~"

I

8.0

+ a,a~ctionpethod (15 v e c t o r s )

16.0

24.0

I

i [

32.0

UL ET h 4

Fig. 7. Accuracy of postbuckling solutions obtained using multiple-parameterreduction method. Sixteen-ply graphite-epoxy plate shown in Fig. 4. Roman numerals indicate points of generating basis vectors, q,. = ( N~,L "-/ET,h3).

third and fourth-order derivatives of {X} with respect to A~ and A2) are less than 0.05% and the errors in the second eigenvalue are less than 0.20%. In order to obtain a nonlinear solution in the vicinity of the stability boundary, the buckling mode shape for a compressive loading equal to 0.715 of the buckling load (N~ = 0.715 N~¢,) was used as a predictor and a corrected solution was obtained using Newton-Raphson iterative technique at q2=(N~2L2/Erha)=26.863. A total of fourteen path derivatives were also generated (these include all the first, second, third and fourth-order derivatives of {X} with respect to q~ and q2). The vectors were orthonormalized using the Gram-Schmidt technique in order to improve the conditioning of the resulting reduced equations. The fifteen basis vectors (the nonlinear solution and the fourteen path derivatives) were used to trace the postbuckling paths corresponding to a fixed compressive loading and an increasing shear load. The postbuckling curves for the three cases of N~ =0.25, 0.50 and 0.715N~cr are shown in Fig. 7. The fifteen vectors were used to advance the solution corresponding to N~ =0.7!5N~c, until q: =77452 corresponding to wah = 2.9. The maximum principal strain at that loading was 5.61 × 10-'. The same set of basis vectors was used to trace the postbuckling paths for the other two cases corresponding to N~ = 0.25 and 0.5ON,, until wc/h = 3.0 (see Fig. 7). As can be seen from Fig. 7, the transverse displacements and strain energies predicted by the reduced system are highly accurate up to w<,/h = 2.0 (corresponding to a maximum principal strain of 3.17× l0 -3 for the case of N~ =0.715N~c,). As wc/h increases beyond 2.0, the predictions of the reduced system start to deviate from those of the full system. The maximum errors in the transverse displacement w¢ obtained by the reduction method were 5.51% for N~ = 0.25N~,, q2 = 86.57; 1.16% for N~ = 0.5N~,. q-_= 83.15, and 4.22% for N~=0.715N~,, q2=77.45. The corresponding errors in the total strain energy (measured at the same values of q, and q~) were 7.28, 8.70 and 8.48%.

76

~.

K. Noor and ]. M. PfTErS

ix, 2

E = 6,895 = 101ON/re2

A = 6,5 ~ 10 "4 m 2

~xl

L = 6.0 in

f = 0.60 m

{ : 0.;'5 m

Equation of Surface 2

t

t

+ x2 +

2

+ 7.2) 2 = 6 0 . 8 4

Loading, Systems

x3

ql;

f

= Xl i

(x3

-'L

concentrated loads in x 3 direction at all nodes concentrated loads in ~3 a t nodes having x 2 < 0

q2:

~i

direction

Fig. 8. Geodesic dome used in the present study.

5.3 Geodesic dome subjected to combined loading The next problem considered is that of determining the stability boundary of a shallow geodesic dome subjected to two independent conservative loading systems. The geometric and material characteristics of the dome is shown in Fig. 8. The loading systems consist of equal concentrated loads acting in the x3-direction. In the first system the entire dome is loaded and in the second, only the nodes whose coordinate x,. < 0 are loaded. The two loading systems are characterized by the two parameters q, and q2, which are selected to be the magnitudes of the nodal concentrated loads. The dome has 37 internal nodes (a total of 111 nonzero displacement degrees of freedom and 132 nonzero force

.......

I O +

×

-d.2 [

103 q

unknowns). Twenty basis vectors were generated for the unloaded dome (q, = q2=0, {X}=0, and {H}=O), and the reduced equations were generated. The same set of basis vectors was used for determimng the stability boundary and tracing the nonlinear response of the dome associated with various combinations of ql and q=. The high accuracy of the predictions of the reduction method for three different values of ql is demonstrated in Fig. 9. The stability boundary of the dome was obtained by selecting a value of the parameter ql and incrementing the other parameter q2 until a critical point is encountered (corresponding to a change in the sign of the current stiffness parameter S, eqn (18)), then repeating the process using other values of q~. The accuracy of the

Full system

14 vectors i Reduction 20 vectors ~ rneth~ III II 71[

= 2 3 ,< 1037

,x 2

-3.15

J V V ' ~ / x Y V \ yM / 3

QZ

~

~,~

ql = -4.6 × lo ~

Newtons

-2.10

JVVV~V'~Y',,_

.,.,,,...,,~r, ~

'~j~/'~/X/x/'X/X/X~-

~

x '

I

"~J\/~A,/V',A,,E~

-

-t.05

// / 0 3

=-6.9×1 0,~

-0 g

-;.2

,va ,n meters

Fig. 9. Accuracy of solutions obtained by muRiple-parameter reduction method used in conjunction with the mixed formulation. Geodesic dome shown in Fig. 8. q, are concentrated loads at all nodes, q, are concentrated loads at nodes having x: < 0.

77

Recent advances in reduction methods for instability analysis of structures

I "~N

I

Full syitem 1 14 vectors ~ Reduction I

0

,+

20 vectorsI

rneth~

l

x2

j-\,/~ J\,IX,,°\/',,,.

q2

-3.0

~ ~

./ / v ~ / V v v ' ~ ./X/VL/VX,/'~/V\..., "_X/X/X/'X/'J ,~,~,./~ '" -',,,/'X/~/~/~,/X/'as/"

,

~I

"J'v'VX/,,/~/"

Newto.s

~x3 -1.5

" ~ ] ~ ~ -

Xl

\ D

',

',

1

-2.5

-5 0

-7.5

0.7"~

I × 103 -I0

ql in Newtons Fig. 10. Accuracy of stability boundary obtained by multiple-parameter reduction methodused in conjunction with the mixed formulation. Geodesic dome shown in Fig. 8. q, are concentrated loads at all nodes, q: are concentrated loads at nodes having x2 < 0.

~ 6 . 8 9 5 x 1010 N / m 2 0.3

~

R t 2,54

[

*

10 "2 rn

~..~.,,..~.....

h = 2.54 ~' 10 -4 rn x 10 -4 m 2

A z 2.54

h3

It1- ~

(a) Closed rm9

(b) Semicircular arch

Fig. 11. Lowest buckling modes for circular ring and semicircular arch subjected to hydrostatic pressure loading.

Table 2. Accuracy and convergence of minimum eigenvalues ~ = p(R~IEI) obtained using reduction method, Circular ring and semicircular arch subjected to uniform pressure loading (see Fig. 11). Number of Basis Vectors

Circular R/n8 Co~ta~t iDirection

3.9997 3.9997

Hydrostatic

Semicircular Arch Cobstant Direction

I

Bydrostatic

2.9997

3.2720

8.0006

2.9997

3.2709

6.7224

3.2709.

3.0037 3.0000 2.9997

Full System

3.9997

2.9997

3.2709

2.9997

', K Nol)Rarid i ~t. PETERS stability boundary obtained using the multiple-parameter reduction method with 14 and 20 vectors is shown in Fig. 10. As can be seen from this figure, the 14 basis vectors do not accurately predict the critical points associated with large values of q~ ii.e. for q, close to q ~ and q,_ =0). On the other hand. the twenty basis vectors accurately predict the entire stability boundary.

! od.

:4

where d is a typical design variable, rhe ~ensitivlt~ coefficients {OXlOd} are then obtained by differentiating the transformation relations, eqns (4), with respect to d: aF

¸

,

5.4 Bifurcation buckling of circular ring and semicircular arch subjected to hydrostatic pressure The last set of problems considered is that of bifurcation buckling of circular ring and semicircular arch subjected to hydrostatic pressure. The material and geometric properties of both structures are shown in Fig. 1t. Analytical solutions for both problems are given in [45, 46]. Finite element solutions are given in [41,47] for the two cases of constant direction and hydrostatic loadings. Due to symmetry, only one quadrant of the ring was considered and was modeled by using four shear-flexible deep curved-beam elements with quintic interpolation functions for each of the displacement and rotation components (a total of 59 nonzero degrees of freedom). The arch was modeled using eight shear-flexible elements (a total of 119 nonzero degrees of freedom). The basis vectors were generated for the unloaded structures (q = 0, {X} = 0), then the reduced stiffness. geometric stiffness and load stiffness matrices were generated and used to compute the critical values of the hydrostatic pressure p , . The buckling mode shapes for both the ring and the arch are shown in Fig. 11. The convergence of the bifurcation buckling load p¢~ with the increase in the number of basis vectors is shown in Table 2 for the two cases of hydrostatic and constant direction pressures. In the latter case. the load stiffness matrix is neglected. As can be seen from Table 2. the buckling pressure for the ring obtained using three basis vectors agree to five significant digits with the full system solutions. For the semicircular arch similar agreement is obtained by using seven basis vectors.

6. FUTUREDIRgCTIONSFOR RESEARCH IN ~DOCTION Mm"aODS Among the different research areas which have high potential for application of reduction methods are the following:

6.1 Use of reduction methods in sensitivity calculations In automated optimum design applications, sensitivity of the response of the system to changes in the design variables is needed. This information is provided by calculating the derivatives of the response quantities with respect to design variables. These are referred to as sensitivity calculations. For linear systems, several efficient techniques have been proposed for sensitivity calculations (see, e.g. [48-50]). However, only few applications have been reported for nonlinear systems. Reduction methods appear to have high potential for sensitivity calculations. A possible approach for applying these techniques is to differentiate the reduced nonlinear equations, e.g. eqns (6), with respect to the design variables as follows:

The coefficients of the matrix [aF/~d] are obtained by differentiating each of the recurrence relations for obtaining the basis vectors with respect to d. For the effective use of this approach in sensitivity calculations, an efficient computational algorithm needs to be developed to reduce the effort in generating

{aX/ad}. 6.2 Use o]: reduction methods in mode interaction prob-

lems Mode interaction problems have attracted considerable attention in the last decade. A review of the work done on this subject up to 1976 is given in [51]. There are basically two types of mode interaction problems. The first involves the interaction of recognizable local and overall modes. An example of this is the interaction between overall Euler buckling and local member buckling of latticed columns. The second type of modal interaction is related to local imperfections in the structure which have the effect of decreasing the stiffness in such a way as to decrease the critical load corresponding to general instability. An example of the second type of modal interaction is provided by the interaction between the column-type modes and the local modes in stiffened prismatic structures subjected to axial compressive loading. Instability analysis in the presence of mode interactions generally requires the use of a refined model for the initial discretization and is, therefore, quite expensive. Application of reduction methods to these problems can result in substantial savings in the computational effort. L CONCLUDINGI~MARKS A number of aspects of reduction methods for instability analysis of structures are discussed including: (a) multi#e-parameter reduction methods for instability analysis of structures subjected to combined loads; (b) application of reduction methods to instability problems with displacement-dependent and nonconservative toadings; and (c) use of mixed models in conjunction with reduction methods for instability analysis. Four numerical examples are presented to demonstrate the effectiveness of using reduction methods for instability analysis of structures. The four problems are: (a) instability of shallow spherical cap subjected to combim~l external pressure and concentrated load. (b) bifurcation buckling and postbuclding of laminated composite plate subjected to combined in-plane compressive and shear Ioadings, (c) instability of a geodesic dome with two independent loading parameters: and (d) bifurcation buckling of a circular ring and a semicircular arch subjected to hydrostatic pressure loading. Also, a number of research areas which appear to have high potential for application of reduction methods are identified. ]'he results of the present study suggest several conclusions relative to the effectiveness of reduction

Recent advances in reduction methods for instability analysis of structures methods in instability problems and the particular choice of the basis vectors in these problems. (!t The proposed single-parameter and multipleparameter reduction techniques for instability analysis of structures are hybrid methods which combine the major advantages of finite element technique, classical Rayleigh-Ritz technique, and perturbation method, namely: (a) Modeling versatility. (b) Reduction in total number of degrees of freedom. (c) Simplicity of determining the stability boundary and tracing pos1-critical point paths. Moreover, it greatly alleviates the following major drawbacks of the aforementioned three techniques: (a) Excessive amounts of computer time required for the finite element instability analysis of large complex structures. (b) Difficulty of selecting global approximation functions for classical Rayleigh-Ritz technique. (c) Small radius of convergence of the Taylor series expansion used in classical perturbation techniques. (2) The particular choice of the basis vectors in the multiple-parameter reduction method as the derivatives of the nodal displacements (and stress parameters) with respect to the independent path (or control) parameters allows the use of the same set of basis vectors and same system of reduced equations/or: (a) the generation of the entire stability boundar3.,, and (b) the computation of a number of post-bifurcation solution paths (corresponding to different combinations of the load parameters). Typically, in each of the postbuckling paths one of the path parameters is varied and the others are fixed. The generation of diffei'ent postbuckling curves provides a quantitative assessment of the sensitivity of the postbuckling response to changes in the load parameters. (31 The multiple-parameter reduction method can be easily applied to problems with several independent load parameters. However, as the number of independent load parameters increases, both the number of basis vectors that need to be generated and the size of the reduced system of equations increase at a faster rate. Therefore, for practical applications an attempt should be made to minimize the number of independent load parameters by treating some of the load parameters as functions of others. (41 By proper selection of the path (or control) parameters, the proposed multiple-parameter reduction method can be used to study imperfection sensitivity as well as the sensitivity of postbuckling response to changes in the design variables of the structure. (51 The use of reduction methods in conjunction with mixed models for instability analysis of structures offers the following two major advantages over the displacement approach: (al Lower computational expense in generating the basis vectors. (b) Better approximation for the reduced solution. The efficiency of the reduction method can be enhanced by using mixed models with discontinuous stresses at interelement boundaries and eliminating the stress parameters as well as their path derivatives on the element level. 16) Reduction methods are very effective tools for instability analysis of structures subjected to pseudo-

79

static displacement-dependent nonconservative loading. In such case the load stiffness matrices [K'"'] and [/~'tq'] are unsymmetric. However, the generation of the basis vectors involves the decomposition of the symmetric linear global stiffness matrix [K], and due to the small number of the reduced equations, no approximation or sy_mmetrization of the unsymmetric load ~tiffness matrix [K ~"~] is needed. REFERENCES

I. W. T. Koiter, On the stability of elastic equilibrium. Thesis, Delft, H. J. Paris, Amsterdam, i945 (in Dutch), English translations (a) NASA Tr F-10833, 1%7, (b) AFFDL-TR-7025, 1970. 2. J. W. Hutchinson and W. T. Koiter, Postbuckling theory. Appl. Meck Rer. 23, 1353-1366 (1971). 3. R. H. Gallagher, The finite element method in shell stability analysis. Comput. Structures 3, 543-557 (1973). 4. B. Budiansky, Theory of buckling and postbuckling of elastic structures. Advances in Applied Mechanics, VoI. 14, pp. 1-65. (Edited by C. S. Yih). Academic Press, New York (1974). 5. B. O. Almroth and J. H. Starnes, The computer in shell stability analysis. J. Engng Meck Div., ASCE 101(EM6), 873-888 (1975). 6. W. T. Koiter, Current trends in the theory of buckling. Proc. IUTAM Syrup. on Buckling of Structures (Edited by B. Budiansky), pp. 1-16. Springer-Verlag, Berlin (1976). 7. B. Budiansky and J. W. Hutchinson. Buckling: progress and challenge. Trends in Solid Mechanics (Edited by J. F. Besseling and A. M. A. Van der Heijden). Sijthoff-Noordhoff. Amsterdam (1980). 8. D. Bushnell, Buckling of shells--pitfalls for designers. AIAA J. 19, 1183-1226 (1981). 9. A. K. Noor, Survey of computer programs for solution of nonlinear structural and solid mechanics problems. Comput. Structures 13,425.-465 ,~1981). 10. J. Connor and N. Morin, Perturbation technique in the analysis of geometrically nonlinear shells. High Speed Computing of Elastic Structures (Edited by B. Fraeijs de Veubeke), pp. 683-706. University of Liege (1971). ii. R. H. Gallagher, Perturbation procedures in nonlinear finite element structural analysis. Computational Mechanics-Lecture Notes in Mathematics, Vol. 461. Springer-Verlag, Berlin (1975). 12. J. L. Batoz and G. Dhatt, Incremental displacement algorithms for nonlinear problems. Int. J. Num. Meth. Engng 14, 1262-1267 (1979). 13. E. Riks. An incremental approach to the solution of snapping and buckling problems. Int. J. Solids Structures 15. 529-55! (!979). 14. G. P6nisch and H. Schwetlick, Computing turning points of curves implicitly defined by nonlinear equations depending on a parameter. Computing 26, 107-121 (1981~. 15. E. Ramm. Stralegies for tracing nonlinear response near limit points. Nonlinear Finite Element Analysis in Structural Mechanics, Prec. Europe-U.S. Workshop, Ruhr Universitat Bochum, West German~,, 28--31 July, 1980, pp. 63-89. Pergamon Press, Oxford 11981't. 16. M. ~.. Crisfield, A fastiincremental iterative solution procedure that handles "snap-through'. Comput. Structures 13, 55-62 (1981). 17. R. G. Melhem and W. C. Rheinboldt, A comparison of methods for determining turning points of nonlinear equations. University of Pittsburgh, Institute for Computational Mechanics and Applications, Techn. Rep. ICMA 82-32, ]an. 1982. 18. J. E. Dennis and 2. J. More, Quasi-Newton methods, motivation and theory. SIAM Rev. 19.46--89 (1977). 19. H. Matthies and G. Strang, The solution of nonlinear finite element equations. Int. J. Num. Meth. Engng 14. 1613-1626 (1979). 20. K. J. Bathe and S. Ramaswamy, An accelerated subspace

80

21. 22. 23. 24. 25.

26. 27.

28. 29. 30.

31. 32. 33.

34. 35.

~. K, NOORand J, M. PETERS iteration method. Comput. Meth. AppL Mech. Engng 23, 313-331 (1980). A. K. Noor, M. D. Mathers and M. S. Anderson, Exploiting symmetries for efficient postbuckling analysis of composite plates. AIAA J. 15. 24-32 (1976). B. O. Almroth, F. A. Brogan and P. Stern, Automatic choice of global shape functions in structural analysis. AIAA L 16, 525-528 (1978). D. A. Nagy, Modal representation of geometrically nonlinear behavior by the finite element method. Comput. Structures 10, 683-688 (1979). A. K. Noor and J. M. Peters, Reduced basis technique for nonlinear analysis of structures. AIAA J. 15, 455.-462 (1980). A. K. Noor, C. M. Andersen and J. M. Peters, Global-local approach for nonlinear shell analysis. Proc. Seventh ASCE Conf. on Electronic Computation, pp. 63,1-657. Washington University, St. Louis, Missouri, 6-..8 August, 1979. D. A. Nagy and M. K6nig, Geometrically nonlinear finite element behavior using buckling mode superposition. Cornput. Meth. AppL Mech. Engng 19, 447--484 (1979). B. O. Almroth, P. Stehlin and F. A. Brogan, Use of global functions for improvement in efficiency of nonlinear analysis. Proc. AIAA/ ASMEI ASCEI AHS 22rid Structures, Structural Dynamics and Materials Conf., Part 1, pp. 286-292, 6-8 April, 1981. A. K. Noor. Recent advances in reduction methods for nonlinear problems. Comput. Structures 13, 31-4,1 (198t). A. K. Noor and J. M. Peters, Tracing post-limit-point paths with reduced basis technique. Comput. Meth. Appl. Mech. Engng 28. 217-240 (1981). A. K. Noor and J. M. Peters, Bifurcation and postbuckling analysis of laminated composite plates via reduced basis technique. Comput. Meth. AppL Mech. Engng 29, 271-295 (1981). K. Huseyin, Nonlinear Theory. of Elastic Stability. Noordhoff, Leyden, The Netherlands (1975). H. Stumpf, On the linear and nonlinear stability analysis in the theory of thin elastic shells, lngenieur-Archiv 51, 195-213 (1981). A. Maewal and W. Nachbar, Stable postbuckling equilibria of axially compressed, elastic circular cylindrical shells: a finite element analysis and comparison with experiments. I Appl. Mech., Trans. ASME 44, 475-481 (1977). P. G. Bergan, G. Horrigmne, B. Kr:Iketand and T. H. S~reide. Solution techniques for nonlinear finite element problems. Int. J. Num. Meth. Engng 12, 1677-1696 (1978). A. K. Noor and J. M. Peters, Nonlinear analysis via globallocal mixed finite element approach. Int. J. Num. Mete Engng 15, 1363-1380 (1980).

36. A. K. Noor and C. M. Andersen, Mixed modeis arid reduced/selective integration displacement models ~'or nonlinear shell analysis. Nonlinear Finite Element Analysis ,,( Plates and Shells, ASME, AMD 48. 119-t46 1198t'L 37. G. A. Cohen. Conservativeness of a normal pres,,ure fieid acting on a shell. AIAA/. 4, 1886 (1966). 38. H. D. Hibbitt, Some follower forces and load stiffnesses, int. J. Num. Meth. Engng 14, 937-941 ~1979). 39. K. Loganathan, S. C. Chang, R. H. Gallagher and J. ~=. Abel. Finite element representation and pressure stiffness in shell stability analysis. Int. J. Num. Meth. Engng 14, 1413-1420 ¢1979). 40. H. A. Mang, Symmetricability of pressure stiffness matrices for shells with loaded free edges. Int. J~ Num. Meth. Engng 15, 981-990 (1980). 41. H. A. Mang and R. H. Gallagher, Finite element analysis of thin shells of general form for displacement-dependent loads. Nonlinear Finite Element Analysis of Plates and Shells, ASME, AMD 48, 65-82 (1981). 42. V. V. Bolotin, Nonconservative .Problems of the Theory of Elastic Stability (in Russian), 1961 translated into English and published by Pergamon Press, Oxford (1965). 43. M. J. Sewell, On configuration-dependent loading. Arch. Rational Mech. Anal. 23, 327-351 (1%7). ,.t4. T. C. Lop and R. M. Evan-lwanowski, Interaction of critical pressures and critical concentrated loads acting on shallow spherical shells. L AppL Mech., ASME. 33, Series E, 612-616 (1966). 45. S. R. Bodner, On the conservativeness of various distributed force systems. J. Aeronaut. Sci. 25, 132-133 (1958). 46. S. P. Timoshenko and J. M. Gere, Theory o[ Elastic Stability, 2nd Edn. McGraw-Hill. New York (1%1). 47. J. L. Batoz, Curved finite elements and shell theories with particular reference to the buckling of circular arch. Int. J. Num. Meth. Engng 14, 774-779 (1979). 48. A. K. Noor and H. E. Lowder, Approximate reanalysis techniques with substructuring. J. Structural Div., ASCE 101, 1687-1698 (1975). ,t9. E. J. Haug and J. S. Arora, Design sensitivity analysis of elastic mechanical systems. Comput. Meth. AppL .~lech. Engng 15. 35-62 (1978). 50. J. S. Arora and E. J. Haug, Sensitivity analysis and optimization of large scale structures. Control and Dynamic Systems, pp. 247-275. Academic Press, New York (1979). 51. V. Tvergaard, Buckling behavior of plate and shell structures. Theoretical and Applied Mechanics (Edited by W T. Koiter), pp. 233--247. North-Holland, Delft (1976).