Mixed formulation in Koiter analysis of thin-walled beams

Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

www.elsevier.com/locate/cma

Mixed formulation in Koiter analysis of thin-walled beams Giovanni Garcea * Dipartimento di Strutture, Universit a della Calabria, 87030 Arcavacata di Rende, Cosenza, Italy Received 21 April 1999; received in revised form 1 December 1999

Abstract The paper extends the asymptotic mixed formulation proposed in (G. Garcea, G. Salerno, Sanitizing of locking in Koiter perturbation analysis through mixed formulation, Fourth World Congress on Computational Mechanics, Buenos Aires, Argentina, 1988; G. Garcea, G. Salerno, R. Casciaro, Comput. Methods. Appl. Mech. Engrg. 180 (1±2) (1999) 137±167) to structures composed of an assemblage of ¯at slender elastic panels. The same type of structure has already been investigated in (A.D. Lanzo, G. Garcea, R. Casciaro, Int. J. Numer. Methods Engrg. 38 (1995) 2325±2345; A.D. Lanzo, G. Garcea, Int. J. Numer. Methods Engrg. 39 (1996) 3007±3031) with a compatible formulation and using the usual hypothesis of negligible precritical rotations in order to avoid the extrapolation locking (see also (G. Garcea, G.A. Trun®o, R. Casciaro, Comput. Methods. Appl. Mech. Engrg. 165 (1±4) (1998) 137±167) produced by the high axial to ¯exural sti€ness ratio of the panels. The use of a mixed formulation provides a consistent way of eliminating the extrapolation locking without introducing ad hoc approximations. It is thus both reliable and accurate. When combined with a quadratic strain±displacement relationship, as presented here, it gives zero fourth-order energy variations, further computational simpli®cations and a `natural' linearization of the eigenvalue problem required by the buckling analysis. An extensive numerical testing shows the e€ectiveness and robustness of the proposed approach. Ó 2001 Elsevier Science B.V. All rights reserved. Keywords: Nonlinear elasticity; Asymptotic methods; Multi±modal buckling; Thin-walled beams; Mixed ®nite elements; Accuracy, reliability

1. Introduction The availability of high-strength materials and the need to minimize the weight of structural elements naturally leads to the use of slender structures. For these reasons, a great amount of research has been done since the 60s with regard to both the theoretical and FEM implementation aspects of the analysis of such structures. In particular, great interest has been devoted to structures composed of an assemblage of slender panels (box-girders, reinforced panels, thin-walled structures), due to their technical relevance for the aeronautical and civil industry. Such structures, due to the high optimization of the design which tends to produce nearly coincident buckling loads, are characterized by coupled instability phenomena and strong imperfection sensitivity. This noticeably complicates the analysis, which has to take into account both the degradation of limit-loads due to additional imperfections and the possible residual strength in the post-buckling range. Thus powerful, that is, reliable, accurate, robust and fast, analysis algorithms are required. The standard industrial approach, based on the so-called path-following method [4,5], while potentially robust and highly accurate, does not appear to be completely satisfactory. In fact, being designed to

*

Tel.: +39-0984-494033; fax: +39-0984-49-4045. E-mail address: [email protected] (G. Garcea).

0045-7825/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 0 ) 0 0 2 6 8 - 1

3370

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

reproduce a single loading process, it requires a large number of successive runs for evaluating the reduction, in load capability, due to all the possible imperfections. Due to the high cost of each run, this usually prevents any optimal design process based on a global safety analysis. Furthermore, we have to consider that, in the presence of multiple, nearly coincident, critical modes, the path-following strategy can become unreliable due to the diculty of recognizing and exactly reconstructing successive bifurcations and modal jumping phenomena which usually occur in such cases. Koiter's asymptotic approach [1±3] represents a suitable alternative, being theoretically able to provide, even in cases of multiple, nearly coincident, buckling modes a simple and complete description of the structural behavior emerging from the bifurcation clusters, including a priori information on the most dangerous shapes of the external imperfection, as described by the so-called `attractive paths theory' [6±12]. In some recent papers [13±16] it has been shown that it is possible to derive from Koiter's theory a ®nite element implementation, suitable for reconstructing structural behavior presenting both bifurcation or snapping. The algorithm is based on the following steps: (i) the fundamental path is obtained as a linear extrapolation starting from the initial known rest con®guration; (ii) a bifurcation search is performed along the fundamental path obtaining the relevant buckling loads and associated modes; (iii) an extrapolation of the total potential energy starting from the ®rst bifurcation point (or from an average point in the bifurcation cluster) is obtained through a consistent asymptotic expansion in terms of modal amplitudes ni ; (iv) the equilibrium path is then obtained by solving a nonlinear system of equations in terms of ni using a standard path-following method. The fundamental idea of the perturbative method is to transform a snapping problem into a bifurcation one, to solve the latter and, ®nally, to reconstruct in an asymptotic fashion, the solution of the snapping problem [13,14,22]. The proposed asymptotic method then furnishes the equilibrium path of a structure for a given loading process, at a fraction of the cost of that provided by the path-following method (essentially that of a linear stability analysis). In addition, once the ®rst analysis has been carried out, successive re-analysis taking into account di€erent imperfections in the geometry or in the loads can be performed with a negligible extra cost. This renders the method suitable for safety analysis based on probabilistic procedures [12]. Note that the description of the nonlinear response of the structure is reduced to a limited number of equations that represent the most meaningful part of the structural behavior (the method achieves the same goal as the socalled modal reduction method [35] but the higher articulation of the projection manifold provides more accurate results). This allows the use of sophisticated path-following strategies in step (iv) and avoids diculties related to secondary bifurcations and modal jumping (e.g. see [10]). A careful implementation of the proposed asymptotic analysis can actually furnish highly-accurate results, as has been shown by both theoretical investigations [22] and numerical tests [14±16]. As the proposed asymptotic formulation is based on an extrapolation, it is particularly sensitive to the representation used and, in particular, on which variables (e.g. displacement and/or stress components) are used for describing the problem and then, directly subjected to the extrapolation process (we call them `primary variables' in [23,25]) and which variables are derived from the ®rst ones using the internal relationship (e.g. strains and stresses obtained from the extrapolated displacements using compatibility equations and the constitutive elastic relationship). The equations at hand being nonlinear, di€erent choices, while equivalent from a strictly asymptotic point of view, can produce very di€erent results at a ®nite distance from the extrapolation point. In a previous paper [25] we showed that unsuitable choices for the primary variables (e.g. those used in a purely compatible formulation in displacements alone) represent the extrapolated equilibrium equation in a highly nonlinear form. The consequent bad representation of the extrapolated fundamental path produces a pathological increase in the values of the bifurcation points and a deterioration in the evaluation of some of the higher-order energy terms which de®ne the post-critical behavior. As a result, the accuracy of the equilibrium paths provided by the method is completely destroyed. In [25] we called this pathological phenomenon `extrapolation locking' and showed that it is essentially due to an interaction between the strong sti€ness ratio (extensional/¯exural) of the elements and the presence of even

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

3371

small rigid rotations of the elements along the fundamental path. We also showed two possible remedies for the problem: (i) To ignore the rotations produced along the fundamental path, i.e. to assume that the bifurcation geometry coincide with the initial one. This hypothesis, usually called `frozen con®guration', is frequently used in buckling analysis and represents an acceptable approximation for a large class of problems [26±31]. (ii) To use a mixed formulation, i.e. to assume as primary variables both displacements and stresses. This decouples the evaluation of the stresses from the rotations and, as a consequence, allows the locking phenomenon to be overcome without introducing additional approximations (see also the discussion given in [23] in a slightly di€erent context). The ®rst method was used in some previous papers [15,16]. Good results in the evaluation of the equilibrium path for a large number of problems were obtained, in excellent agreement with `reference' solutions obtained by carefully tuned path-following analysis. However, in some cases characterized by a quite strong nonlinearity of the precritical behavior, an appreciable di€erence was experienced between the reference results and those provided by the asymptotic method. This suggested a re®nement of the analysis by introducing a mixed formulation. This paper reports the results of this research showing that it is possible to avoid extrapolation locking and also obtain some additional computational advantages. The proposed algorithm proves to be extremely robust and suitable for accurately evaluating equilibrium paths even in problems characterized by a strongly nonlinear precritical behavior. The paper is organized as follows: Section 2 describes the proposed mixed formulation for plate assemblages; Section 3 presents a brief sketch and the principal formula of the asymptotic method used; Section 4 describes the proposed ®nite element discretization, including all the relevant implementation details; Section 5 discusses the results obtained in extensive numerical testing; some remarks and ®nal comments are then given in the concluding section.

2. The mixed formulation for plate assemblages As in [15,16] we consider thin-walled structures that can be represented as assemblages of rectangular plates interconnected along their longitudinal edges (see Fig. 1), under general boundary conditions, and subjected to a load increasing linearly according to an ampli®er k. Each plate is assumed to be constant in thickness, which is small compared with the length and the width of the plate. For each plate a local Cartesian coordinate system …x; y; z† is assumed, where the x±y plane describes the undeformed middle

Fig. 1. Typical cross-sections with related vertexes and panels.

3372

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

surface. In such a system …u; v; w† will be the components of the displacement ®eld u‰x; yŠ of the middle surfaces measured from the undeformed con®guration. 2.1. The plate model In previous papers, see [20,21], it has been shown that the asymptotic approach requires careful kinematical modeling in order to obtain a consistent evaluation of the third- and fourth-order energy terms used in the analysis. This implies that the kinematical modeling must be able to exactly evaluate the fourth-order derivatives with respect to the displacements. In general this is not obtained with technical theories and in particular with the technical plate models which, on the base of a quadratic strain±displacements relationship, zero all the third and fourth-order displacement terms usually smaller of the second-order ones. This produces a kinematical inconsistency in the evaluation of the postcritical curvature de®ned as a di€erence between fourth-order and second-order energy variations. In a plane frame context [20], it has been shown that the contribution given by the second-order derivatives of the strain tends to vanish and then the exact evaluation of the overall behavior is de®ned only by the high order kinematical terms. In [15,16], we showed that the classical Karman±Marguerre nonlinear plate theory, enriched with the appropriate nonlinear terms related to the in-plane rotations, is suciently accurate, for the structures we want to analyse due to the strong stress redistribution which follows the buckling. This redistribution makes the quadratic terms, in the energy variation, di€erent from zero and notably bigger than the omitted high order ones. This e€ect, typical of three-dimensional structures built up with several plates [15,16,19,20,33], makes the error regarding an exact evaluation of the postcritical curvature negligible. However for thin-walled beams of oblong form, i.e with a high ratio between the beam axis length and the maximal cross-section dimension, presenting global displacements (i.e. minimal local distortion) this three-dimensional constraint e€ect can be completely eliminated. For these structures the use of proper kinematics [20,33] is a fundamental prerequisite to obtain an accurate description. In [16] we have investigated two technical models in particular. The ®rst, that we called Simpli®ed Lagrangian, is characterized by the following strain de®nition: 1 ex ‰uŠ :ˆ u;x ‡ …v2 ;x ‡w2 ;x †; vx ‰uŠ ˆ w;xx ; 2 1 2 ey ‰uŠ :ˆ v;y ‡ …u ;y ‡w2 ;y †; vy ‰uŠ ˆ w;yy ; 2 cxy ‰uŠ :ˆ u;y ‡v;x ‡w;x w;y ; vxy ‰uŠ ˆ w;xy :

…1†

The second, called Compete Lagrangian, only differs in the in-plane strain relationship, de®ned as follows: 1 ex ‰uŠ :ˆ u;x ‡ …u2 ;x ‡v2 ;x ‡w2 ;x †; 2 1 2 ey ‰uŠ :ˆ v;y ‡ …u ;y ‡v2 ;y ‡w2 ;y †; 2 cxy ‰uŠ :ˆ u;y ‡v;x ‡u;x u;y ‡v;x v;y ‡w;x w;y :

…2†

In the previous equations the vector u 8 9 < u‰x; yŠ = u :ˆ v‰x; yŠ : ; w‰x; yŠ collects the components in the x, y and z directions of the displacement ®eld, and a comma followed by x; y denotes partial differentiation with respect to x; y. We have also shown that the two kinematical models are in good agreement and describe practically the same behavior when a relevant stress redistribution effect is present as in the present context. It was also shown that there is a small error in both in-plane and out-ofplane plate behavior with the Simpli®ed Lagrangian model. The Complete Lagrangian is accurate when the

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

3373

structure behavior is characterized only by in plane displacements but there is a greater error with regard to the out-of-plane behavior. This leads to the use of the Simpli®ed Lagrangian plate model. In vector form the strain 8 9 8 9 < ex ‰uŠ = < vx ‰uŠ = ey ‰uŠ ; v‰uŠ ˆ vy ‰uŠ ‰uŠ ˆ …3† : c ‰uŠ ; : v ‰uŠ ; xy xy and the stress resultants 8 9 8 9 < Nx ‰x; yŠ = < Mx ‰x; yŠ = n ˆ Ny ‰x; yŠ ; m ˆ My ‰x; yŠ : : ; ; Nxy ‰x; yŠ Mxy ‰x; yŠ

…4†

are related to each other through the following linear constitutive relations: n ˆ C;

m ˆ Dv;

 ˆ Fn:

…5†

The elastic 3  3 matrices C and D (F ˆ C 1 ) are symmetric and positive de®nite. By indicating with u the con®guration of the structure, the strain energy of the plate is expressed, in compatible form, by Z n o 1 U‰uŠ :ˆ ‰uŠT C‰uŠ ‡ v‰uŠT Dv‰uŠ dA …6† 2 A or, in mixed form, by Z n 1 U‰uŠ :ˆ nT …‰uŠ 2 A

o T T Fn† ‡ ‰uŠ n ‡ v‰uŠ Dv‰uŠ dA:

…7†

Note that Eq. (6) de®nes the energy as a functional of the displacement ®eld u‰x; yŠ, that is, the con®guration u refers to structural displacements alone, n‰x; yŠ and m‰x; yŠ being de®ned through Eqs. (5) and (7) de®nes the energy as a function of u‰x; yŠ and n‰x; yŠ, and in this case u refers to both the displacements and the in-plane stresses. 2.2. Relevant variations of the strain energy Denoting the Frechet derivation with respect to the con®guration u by a prime, the successive variations of the compatible form of the energy are Z n o 00 _ ˆ _ T Cd‰uŠ ‡ v_ T Ddv dA; U ‰uŠudu nT d_ ‡ ‰uŠ …8† A

Z n o ^‰uŠT Cd_ ‡ ‰uŠ _ ˆ _ T Cd^ ‡ d‰uŠT C^_ dA; U000 ‰uŠ^ uudu A

_ ˆ U0000 u~u^udu

Z n o ~ ^T Cd_ ‡ ~_ T Cd^ ‡ d~T C^_ dA: A

For the mixed form of the energy these variations are: Z n o 00 _ ˆ _ ‡ v_ T Ddv dA; U ‰uŠudu nT d_ ‡ n_ T …d‰uŠ Fdn† ‡ dnT ‰uŠ A

_ ˆ U000 u^udu

Z n o ^ nT d_ ‡ n_ T d^ ‡ dnT ^_ dA; A

_ ˆ 0: U0000 u~u^udu

…9† …10†

…11† …12† …13†

3374

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

_ d, dv, d, _ v, _ n_ and dn are de®ned by: In both cases, the quantities , 8 9 8 9 _ xx = < w; < dw;xx = _ yy ; dv :ˆ v0 ‰uŠdu ˆ dw;yy ; v_ :ˆ v0 ‰uŠu_ ˆ w; : ; : ; _ xy dw;xy w;

…14†

8 <

9 _ x ‡v;x v_ ;x ‡w;x w; _ x = u; _ y ‡w;y w; _ y v_ ;y ‡u;y u; ; _ :ˆ 0 ‰uŠu_ ˆ : ; _ y ‡_v;x ‡w;x w; _ y ‡w;y w; _ x u; 8 <

9 du;x ‡v;x dv;x ‡w;x dw;x = dv;y ‡u;y du;y ‡w;y dw;y d :ˆ 0 ‰uŠdu ˆ ; : ; du;y ‡dv;x ‡w;x dw;y ‡w;y dw;x

…15†

8 9 _ x = < dv;x v_ ;x ‡dw;x w; _ y ‡dw;y w; _ y ; _ ˆ du;y u; d_ ˆ 00 ‰uŠudu : ; _ y ‡dw;y w; _ x dw;x w;

…16†

8 9 < N_ x = n_ ˆ N_ y ; : _ ; Nxy

…17†

8 9 < dNx = dn ˆ dNy : : ; dNxy

Note that, in the mixed formulation, the in-plane stresses n correspond to the primary parameters of the con®guration u. In contrast to the compatible formulation the in-plane stress ®eld is derived as a function of the displacement u through the constitutive law n ˆ C. This leads, due to the quadratic deformation± displacement relationship used, to an expression of the strain energy that is cubic in u for the mixed approach. Consequently the third variation is independent of the con®guration and the fourth variation is zero while for the compatible formulation the third variation is linear in u and so the fourth variation is constant and in general other than zero. The greater simplicity of the mixed formulation with its lower order energy representation is evident. 3. The asymptotic method In this section we recall the asymptotic method as proposed in [14±21], for solving both bifurcation and snapping problems, by reporting its relevant equations and procedures.We consider a slender hyperelastic structure subjected to conservative loads k^ p increasing with the proportionality factor k. The equilibrium is expressed by the virtual work equation U0 ‰uŠdu

^ ˆ 0 8du 2 T; kpdu

…18†

where u 2 U is the ®eld of con®guration variables, U‰uŠ denotes the strain energy, T is the tangent space of U at u and a prime is used for expressing the Frechet derivative with respect to u. Eq. (18) de®nes a curve in the …u; k† space, the equilibrium path of the structure, that can be composed of several branches. We are usually interested in the branch starting from an initial known equilibrium point fu0 ; k0 g and without any loss of generality we can consider u0 ˆ 0, k0 ˆ 0. The asymptotic method provides an approximation for the equilibrium path by the following steps: (1) The fundamental path is obtained as a linear extrapolation, from a known equilibrium con®guration uf ‰kŠ :ˆ k^ u0 ;

…19†

where u^0 is the tangent evaluated at fu0 ; k0 g which is obtained as a solution of the linear equation U000 u^du ˆ p^du

8du 2 T …U000  U00 ‰u0 Š†:

…20†

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

3375

(2) A cluster of buckling loads fk1 ; . . . ; km g and associated buckling modes …_v1 ; . . . ; v_ m † are de®ned along this path by the bifurcation condition U00 ‰kŠ_vdu ˆ 0

8du …U00 ‰kŠ  U00 ‰uf ‰kŠ†:

…21†

The buckling loads are considered to be suciently close to each other to allow the following linearization U00b v_ i du ‡ …ki

^v_ i du ˆ 0 kb †U000 bu

8du …U00b  U00 ‰uf ‰kb ŠŠ†:

…22†

kb being an appropiate reference value of k (the ®rst of ki or their mean value). Normalizing we obtain ^v_ i v_ j ˆ U000 bu

dij

…dij : Kronekers0 symbol†:

…23†

Eq. (21), de®nes a nonlinear eigenvalue problem. Robust iterative algorithms are available for nonlinear eigenvalue analysis [32,21], although they tend to slow down in cases of multimodal search. (3) The tangent space T is decomposed into two orthogonal subspaces: ( P V :ˆ f_v ˆ miˆ1 ni v_ i g …critical subspace†; T ˆ V  W; …24† ^v_ w ˆ 0g …orthogonal subspace†: W :ˆ fw 2 T : U000 bu The asymptotic approximation for the required path is de®ned by the expansion u0 ‡ v_ ‡ w‰k; nk Š; v_ 2 V; w 2 W; u‰k; nk Š :ˆ k^ m 1 1X w‰k; nk Š :ˆ k2 w00 ‡ n n wij ; 2 2 i;jˆ1 i j

…25†

where wij are quadratic corrections introduced to satisfy the projection of Eq. (18) onto W and obtained by the linear orthogonal equations U00b wij dw ˆ

_ i v_ j dw; U000 bv

wij ; dw 2 W;

i; j ˆ 0; . . . ; m;

…26†

where v_ 0 :ˆ u^ and w0i ; i 6ˆ 0 are zero due to the orthogonality condition used. (4) The following energy terms are computed for i; j; k ˆ 1; . . . ; m: 1 2 000 2 1 k Ub u^ v_ k ‡ k2 …k 2 6 _ _ _ ˆ U000 ; v v v b i j k

lk ‰kŠ ˆ Aijk

Bijhk ˆ B00ik ˆ B0ijk ˆ Cik ˆ

_ i v_ j v_ h v_ k U0000 b v 0000 2 Ub u^ v_ i v_ k ^v_ i v_ j v_ k ; U0000 b u 00  ik ; Ub w00 w

^3 v_ k ; 3kb †U0000 b u

…27† …28†

 hk U00b … wij w 00  ik ; Ub w00 w

 ih w  ik w  jk ‡ w  jh †; ‡w …29† …30†

where lk are called the implicit imperfection factors, which correspond to a consistent fourth-order expansion of the unbalanced work on the fundamental path (i.e. lk ‰kŠ :ˆ …k^ p U0 ‰k^ uŠ†_vk ). (5) The required equilibrium path is then obtained by satisfying the projection of the equilibrium Eq. (18) onto V. In line with the extrapolations made above we have the following expansion into Taylor series lk ‰kŠ ‡ …kk 1 ‡ …k 2

k†nk

kb …k

kb =2†

m X iˆ1

kb †

m X i;jˆ1

ni nj B0ijk ‡

ni Cik ‡

m 1X 1 n n Aijk ‡ …k 2 i;jˆ1 i j 2

m 1 X ni nj nh Bijhk ˆ 0; 6 i;j;hˆ1

2

kb †

k ˆ 1; . . . ; m:

m X

ni B00ik

iˆ1

…31†

Eq. (31) corresponds to a highly nonlinear system in the m unknowns ni and can be solved using a standard Riks method. The method then reduces the complexity of the nonlinear structural behavior, consistent with the asymptotic framework of analysis, to a very limited number of equations that represent

3376

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

the most meaningful part of the structural behavior. From this point of view it resembles the Noor modal reduction methods [35]. It provides the complete post-buckling behavior of the structure including modal interactions and jumping-after-bifurcation phenomena. Once the ®rst analysis has been performed (steps 1±4), the presence of small additional imperfections expressed by a load eq q~‰kŠ and/or an initial displacement eu u~ can be taken into account in the postprocessing phase (solution of Eq. (31)), with a negligible computational extra-cost, by simply adding the additional terms: lqk ‰kŠ :ˆ

e~ q‰kŠ_vk ;

luk ‰kŠ :ˆ kU000 u^u~v_ k

…32†

to (31). That is, the imperfection term lk ‰kŠ is modi®ed to 1 1 ^2 v_ k ‡ k2 …k lk ‰kŠ ˆ k2 U000 bu 2 6

^3 v_ k ‡ lqk ‰kŠ ‡ luk ‰kŠ: 3kb †U0000 b u

…33†

This allows a full imperfection sensitivity analysis to be performed. Other considerations about imperfection sensitivity analysis can be found in the attractive paths theory discussed in [6±9] and in [10]. It provides a conceptual framework for selecting a narrower range of meaningful imperfections by allowing a full (stochastic) imperfection sensitivity analysis to be performed [12]. 3.1. An extrapolation locking phenomenon Even if the described algorithm has been extensively tested and shown capable of providing excellent results (the theoretical analysis carried out in [22] shows a great accuracy in the evaluation of the limit-load and of the initial postcritical behavior), its implementation in the solution of problems characterized by a nonnegligible nonlinear precritical behavior can, sometimes, furnish unexpectedly bad results. This aspect necessitates further consideration and some preliminary remarks. When we substitute a curve with its Taylor expansion of order n, the obtained extrapolation is affected by an in®nitesimal error O…nn †, if n is the distance from the extrapolation point. However, we do not get any a priori quantitative information on the extrapolation error at a ®nite distance. On the other hand the curve subject to extrapolation is one particular representation of the equilibrium path, different representations of the same path being possible. Extrapolations of the same order based on different representations are equivalent from an asymptotic point of view, but can imply different errors at a ®nite distance. This means that the choice of what representation, that is, the selection of what we called primary variables, i.e. the variables describing the path and directly subject to extrapolation, plays a central role in reducing the error at a ®nite distance. The described asymptotic method is based on two independent extrapolations, those of the fundamental and bifurcated paths. While not obvious, it is the ®rst one which requires the maximum care. In fact, since the rest con®guration …u0 ; k ˆ 0† is an equilibrium con®guration available at the beginning of the analysis, it is the most convenient and natural con®guration to be used as the extrapolation point for the fundamental path. On the other hand, the algorithm requires that the bifurcation con®guration …ub ; kb † be known before the energy terms described in step (5) are calculated. But rest and bifurcation con®gurations, both belonging to the fundamental path, are at a ®nite distance each from the other. Therefore a primary objective to guarantee the reliability of the algorithm is to increase the accuracy of the extrapolation of the fundamental path. This objective can be obtained by selecting the smoothest representation among all the possible representations of the equilibrium path, which corresponds to choosing the most suitable extrapolation variables. To give an idea of how important the representation choice in controlling the extrapolation error, it is sucient to bear in mind the example of a circle, whose linear extrapolation in terms of Cartesian coordinates is a€ected by a strong extrapolation error which completely disappears, when using polar coordinates. As shown in [25], a purely compatible formulation (that is, a formulation where the path is described through the extrapolation of displacements and the stresses are derived from the displacements) can present such a high degree of nonlinearity that renders completely unreliable, at ®nite distance, the asymptotic

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

3377

evaluation of the relevant quantities (second, third and fourth energy variations) required by the analysis. This produces a pathological increase in both the buckling load and the post-critical path curvature, that can destroy the accuracy of the solution. Consequently, very large errors can result, usually in the form of a strong locking e€ect. We called this phenomenon extrapolation locking. In [25] (see also [23]), it was shown that the phenomenon is produced by the interaction between even small precritical rotations and high axial/¯exural sti€ness ratios of the elements. Slender plates are usually characterized by a high axial to ¯exural sti€ness ratio and so are particularly sensitive to this problem. Therefore a reliable solution process requires a strategy to avoid the locking. It is worth mentioning that this form of locking is also produced within path-following analysis where, however, it only implies a deterioration of the convergence of the iterative solution process [23]. As shown in [23±25], to which the reader is addressed for a deeper insight into this subject, a mixed formulation, where the in-plane stresses are also assumed as primary variables, gives a much less nonlinear path representation, allowing a sanitization of the locking problems. 3.2. The frozen geometry hypothesis A cure for the locking phenomenon is obtained by introducing a `negligible rotations' hypothesis along the fundamental path, that is, neglecting all the precritical rotations (and then zeroing all the nonlinear precritical e€ects) in the description of the fundamental path. As displacements are usually very small along the fundamental path, this hypothesis is adequate, but can lead to some loss in accuracy for structures characterized by a nonlinear precritical behavior, as shown in [25]. This hypothesis, which implies a linearized buckling analysis, is frequently used (see [26±31]) and implies that all nonlinear terms deriving from k^ u will be neglected in the energy variations with some simpli®cations in the formula of the asymptotic method [15,16]. In particular we set: w00 ˆ 0;

B0ijk ˆ B00ij ˆ Cij ˆ 0

…34†

so that Eq. (31) is simpli®ed to lk ‰kŠ ‡ …kk

k†nk ‡

m m 1X 1 X ni nj Aijk ‡ n n n Bijhk ˆ 0; 2 i;jˆ1 6 i;j;hˆ1 i j h

k ˆ 1; . . . ; m:

…35†

A further computational advantage of this hypothesis is that we have uŠ ˆ U000 ‡ …k U00 ‰k^

^ k0 †U000 0u

and then the eigenvalue search (21) reduces to a standard linear eigenvalue problem i h ^ v_ du ˆ 0 8du 2 T U000 ‡ …k k0 †U000 0u

…36†

by allowing the use of faster standard solution algorithms. Note that the hypothesis of frozen geometry implies that some terms (related to the displacements on the fundamental path) have to be omitted in the expression of U000 0 and in the high-order energy variations used in Eq. (35) (see [15,16] for details and a further discussion of this). A rough linearization of the eigenvalue problem (21) without taking into account all the related consequences of the frozen geometry can introduce kinematical inconsistencies in the formulation and strongly a€ect the evaluation of high-order energy terms Aijk and Bijhk . This circumstance can destroy the accuracy of the analysis as is discussed in [21] and as is shown by the numerical results presented in Section 5 of this paper. 3.3. Mixed formulation When implemented for the mixed plate formulation described in Section 2, the asymptotic method simpli®es noticeably, due to the disappearance of all fourth-order terms. In fact, the expressions (27)±(29) reduce to

3378

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

1 2 000 2 1 ^2 v_ k ; k U0 u^ v_ k ˆ k2 U000 bu 2 2 _ i v_ j v_ k ; ˆ U000 bv

lk ‰kŠ ˆ Aijk

Bijhk ˆ B00ik ˆ

 hk ‡ U00b … wij w  ik ˆ U00b w00 w

…37† …38†

 ih w  ik w  jk ‡ w  jh †; w Cik ;

…39†

B0ijk ˆ 0:0:

…40†

and Eq. (31) becomes lk ‰kŠ ‡ …kk

k†nk

m m m 1 2X 1X 1 X k ni Cik ‡ ni nj Aijk ‡ ni nj nh Bijhk ˆ 0; 2 iˆ1 2 i;jˆ1 6 i;j;hˆ1

k ˆ 1; . . . ; m:

…41†

Besides this formal simpli®cation, the mixed formulation is advantageous in terms of computational costs and accuracy. First of all, due to the absence of the fourth-order terms, the eigenvalue search de®ned by Eq. (21) corresponds exactly to a linear eigenvalue problem. As a consequence, Eq. (24) holds exactly for ®nite distances between the ki so that there is no need to assume they are close. This quickens the solution process in cases of multimodal search and improves the accuracy of the analysis by eliminating the possible rise in kinematical inconsistencies. Even more important is that Bijhk is now expressed by fourth-order term functions of w alone instead of being determined by a small difference between two large contributions, U0000 v_ 4 and 3U00 w2 , as is required by the compatible approach. Thus the evaluation of this quantity is naturally free from interpolation locking due to an inconsistent matching between in-plane and out-of-plane displacements discretization (see [15,33] for an insight into this problem). In this way the ®nite element modeling is simpli®ed and the need for a stress post-processing ®nite element representation is avoided (see [34]). However, the main advantage is that the mixed formulation implies a much less nonlinear representation of the equilibrium path. It naturally allows us both to avoid the extrapolation locking without the need of ad hoc hypotheses and without introducing additional errors in the solution, and to improve the long distance accuracy of the asymptotic solution (refer to [25,23] for a discussion of these topics).

4. Mixed ®nite element formulation for plate assemblies The ®nite element model adopted for general assemblies of plates strictly follows that described in [15,16]. The structure is split into separate panels connected along the longitudinal edges as shown in Fig. 1, each panel being discretized by a rectangular mesh with sides lx and ly element sides. Continuity is assumed at the inter-element boundaries for both displacements (u; v; w) and their derivatives. The code organization is based on substructures, the single panel being treated as a macro-element of the overall structure. All the energy terms needed by the analysis are computed from the elements through a 2  2 Gauss points integration scheme that provides better eciency in the evaluation of the high-order energy terms and allows the ®lter of the interpolation locking (within the compatible formulation) in the estimate of the postbuckling curvature terms Bijhk . The compatible approach with frozen geometry hypothesis and the mixed formulation have both been implemented in an updated release of the KASP code [37]. 4.1. Displacement discretization The displacement discretization is based on the rectangular High Continuity ®nite element proposed in [36] and already used by the authors in [15,16] in the same structural context. Each component of the displacement is interpolated through quadratic splines as shown in Figs. 2 and 3. u is given on the element by

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

3379

Fig. 2. HC interpolation.

Fig. 3. HC shape function.

u‰n1 ; n2 Š ˆ

3 X

hi ‰n1 Šhj ‰n2 Šuij ;

n1 ˆ

i;jˆ1

x

x `x

;

n2 ˆ

y

y `y

 ;

…n1 ; n2 † 2

 1 1 ; ; 2 2

…42†

where h1 ‰nŠ ˆ

1 8

1 1 n ‡ n2 ; 2 2

h2 ‰nŠ ˆ

3 4

n2 ;

h3 ‰nŠ ˆ

1 1 1 ‡ n ‡ n2 8 2 2

x, y are the coordinates of the the center of the element, `x , `y the element dimensions, and the displacement parameters uij are de®ned, according to Fig. 4, on nodes located in the center of the element. The same holds for v and w. The forced inter-element continuity of the displacements and their derivatives reduces the parameters for this local displacement representation, i.e. we obtain a displacement approximation with C 1 continuity in the domain of the panel using a minimal number of interpolation parameters (approximately one parameter per element for each displacement component). The continuity between panels is performed by means of the map

Fig. 4. Local mesh.

3380

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

u ˆ

ua ‡ ub ; 2

v ˆ

va ‡ vb ; 2

ˆ w

wa ‡ wb ; 2

uˆ

wb

wa Dl

 u) suitable that transforms the variables placed near the edges (see Fig. 5) into boundary variables ( u; v; w; for forcing the C 1 continuity of displacement on the panels interfaces (see Fig. 6). 4.2. Numerical integration A 2  2 Gauss scheme has been used for integrating all the energy terms required by the analysis. Let …xg ; y g † be the position of the Gauss point g within the element e and Ag its weight, the components of the displacement gradient are evaluated in matrix form as u;x ‰xg ; y g Š :ˆ

3 X i;jˆ1

hi ‰ng1 Š;x hj ‰ng2 Šuij ˆ

9 X kˆ1

g Hx;k uk ˆ Hgx ue ;

where the nodal values u‰x; yŠ pertaining to the element are renumbered as k ˆ 3j ‡ i the element vector uTe ˆ ‰u1 ; u2 ; . . . ; u9 Š (see Fig. (4)), and Hgx :ˆ ‰h1 ‰ng1 Š;x h1 ‰ng2 Š; h2 ‰ng1 Š;x h1 ‰ng2 Š; . . . ; h3 ‰ng1 Š;x h3 ‰ng2 ŠŠ

Fig. 5. Simple plate mesh.

Fig. 6. Kinematical relationships along common panel boundaries.

1 and collected in

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

3381

is a ‰1  9Š matrix associated to the Gauss point. In the same manner we de®ne: u;gy ˆ Hgy ue ;

u;gxx ˆ Hgxx ue ; . . . ; w;gyy ˆ Hgyy we :

…43†

Then the element contribution to the quadratic terms appearing in the second variations of the strain energy can be evaluated as Z X X u;x u;x dA ˆ u;gx u;gx Ag ˆ Ag …Hgx ue †…Hgx ue †: …44† Ae

g

g

The contribution to the cubic and fourth-order terms appearing in the third and fourth variations is obtained in the form: Z X u;x w;x w;x dA ˆ Ag …Hgx ue †…Hgx we †…Hgx we †; …45† Ae

g

Z Ae

w;x w;x w;x w;x dA ˆ

X g

Ag …Hgx we †…Hgx we †…Hgx we †…Hgx we †:

…46†

It is evident that the evaluation of the higher order terms involves considerable numerical operations. However they are required only for obtaining scalar quantities or, in the worst cases, vector quantities. In this way their evaluation is not critical with respect to the overall computational costs required by the analysis. 4.3. Stress discretization The in-plane stresses n‰x; yŠ required by the mixed formulation are recovered on Gauss points alone, by using a local criterion for de®ning the constitutive relationships: ng ˆ Cg ;

g ˆ Fng :

ng :ˆ fNx ; Ny ; Nxy g and g :ˆ fex ; ey ; cxy g being the local in-plane stress and strain vectors associated to the single Gauss point. This implies a constant piecewise interpolation of the in-plane element stress. This choice, which may not be optimal for accuracy is justi®ed because: (i) it is the simplest from both the computational and implementational point of view; (ii) we obtain an exact correspondence between compatible and mixed ®nite element interpolation which allows us to concentrate on di€erences in the extrapolation treatment of the stress variables without introducing any disturbing discretization improvement. Still, it is worth mentioning that, while being identical from a local point of view, the mixed and compatible models behave di€erently within the extrapolation process required by the asymptotic analysis. They provide at ®nite distance from the extrapolation point a di€erent response, more accurate in the mixed case, as shown by the numerical results in Section 5. 4.4. Matrix expressions for the energy variations Introducing the following de®nitions (see Fig. 4): de :ˆ fue ; ve ; we g the vector collecting the 18 nodal displacements associated to the element; ne :ˆ fn1 ; n2 ; n3 ; n4 g the vector of the Gauss point in-plane stress variables associated to the element; e :ˆ f1 ; 2 ; 3 ; 4 g the vector of the Gauss point in±plane strain variables associated to the element; ve :ˆ fv1 ; v2 ; v3 ; v4 g the vector of the Gauss point out-of-plane strain variables; Ue :ˆ fne ; de g the local vector collecting all variables related to the element; Ug :ˆ fng ; de g the local vector collecting all variables related to the Gauss point g;

3382

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

we have:  
1 g ˆ Lg ‡ Qg ‰de Š de ; _ g ˆ Lg ‡ Qg ‰de Š d_ e ; 2 d_ g ˆ Qg ‰dde Šd_ e ˆ Qg ‰d_ e Šdde ; v_ g ˆ Sg d_ e ;

…47†

where the following notation is adopted (index g is omitted for easier writing): 2 3 2 3 Hx 0 0 0 0 Hxx Lg ˆ 4 0 Hy 0 5; Sg ˆ 4 0 0 Hyy 5; Hy Hx 0 0 0 Hxy 2

0 6 uT H T H y Qg ‰de Š ˆ 6 4 e y 0

vTe HTx Hx 0 0

wTe HTx Hx wTe HTy Hy

…48†

3

7 7   5: T T T w e Hx Hy ‡ Hy Hx

…49†

The expression of the relevant variations of the strain energy Eqs. (8)±(12) can be written in matrix form, for the compatible formulation _ ˆ U00 ‰uŠudu

o Xn nTg d_ g ‡ _ Tg Cdg ‡ v_ Tg Ddvg Ag e;g

ˆ

X e;g

_ ˆ U000 ‰uŠ^ uudu

n o
T
T ddTe Qg ‰d_ e Š ng ‡ Lg ‡ Qg ‰de Š C Lg ‡ Qg ‰de Š d_ e ‡ STg DSg d_ e Ag ;

…50†

o Xn ^Tg Cd_ g ‡ _ Tg Cd^g ‡ dTg C^_ g Ag e;g

Xn
T
T ^ ˆ dTe Lg ‡ Qg ‰de Š CQg ‰dde Šd_ e ‡ d_ Te Lg ‡ Qg ‰de Š CQg ‰^de Šdde e;g

o
T ‡ ddTe Lg ‡ Qg ‰de Š CQg ‰^ de Šd_ e Ag ; _ udu ˆ U0000 u^u~

o Xn ~ ^g Cd_ g ‡ ~_ Tg Cd^g ‡ d~Tg C^_ g Ag e;g

ˆ

…51†

Xn e;g

o ^ de ŠCQg ‰dde Šd_ e ‡ d_ Te QTg ‰~de ŠCQg ‰^de Šdde ‡ ddTe QTg ‰~de ŠCQg ‰^de Šd_ e Ag dTe QTg ‰~

and, for the mixed one: Xn _ ˆ U00 ‰uŠudu nTg d_ g ‡ n_ Tg dg ‡ dnTg …_ g e;g

…52†

o Fn_ g † ‡ v_ Tg Ddvg Ag

  Xn 
T T ˆ ddTe Qg ‰d_ e Š ng ‡ Lg ‡ Qg ‰de Š n_ g ‡ STg DSg d_ e ‡ dnTg …Lg ‡ Qg ‰de Š†d_ e

Fn_ g

o

Ag ;

e;g

…53† _ ˆ U000 u^udu

Xn e;g

ˆ

Xn e;g

o ^ nTg d_ g ‡ n_ Tg d^g ‡ dnTg ^_ g Ag o ^ nTg Qg ‰dde Šd_ e ‡ n_ Tg Qg ‰^ de Šdde ‡ dnTg Qg ‰^de Šd_ e Ag :

…54†

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

3383

4.4.1. Incremental internal forces vector The element contribution to the incremental vector of internal forces s_ e is de®ned by the energy equivalence: X X _ :ˆ U00 ‰uŠudu dUTe s_ e :ˆ dUTg s_ g : e

e;g

By some manipulations of Eqs. (52) and (54), we can derive, the incremental internal forces vector for the compatible and mixed models: o Xn
T T s_ ce :ˆ STg Dv_ g ‡ Qg ‰d_ e Š Cg ‡ Lg ‡ Qg ‰de Š C_ g Ag ; …55† g

s_ me :ˆ

X

(

g

_ g Fn_ g
T T STg Dv_ g ‡ Qg ‰d_ e Š ng ‡ Lg ‡ Qg ‰de Š n_ g c

) Ag ;

…56†

m

where superscripts …
† and …
† denote the quantities relative to the compatible and mixed formulation. Note that the formal di€erences between the two formulations are limited to the way of calculating the stress quantities ng and n_ g : directly extrapolated for the mixed formulation, derived by the strains through the constitutive equation ncg :ˆ Cg and n_ cg :ˆ Cg for the compatible one. 4.4.2. The tangent sti€ness matrix The element contribution to the tangent sti€ness matrix is de®ned by the following energy equivalence: X X _eˆ _ g: _ :ˆ U00 ‰uŠudu dUTe K_ e U dUTg K_ g U e

e;g

Introducing matrix G‰ng Š 2 0 Ny HTy Hy 6 G‰ng Š ˆ 4 0 Nx HTx Hx 0 0

3 0 7 5 0 Nx HTx Hx ‡ Ny HTy Hy ‡ Nxy …HTx Hy ‡ HTy Hx

de®ned by the following equivalence G‰ng Šd_ e  QTg …d_ e †ng ; we obtain the element sti€ness matrix that for the compatible model can be written as o X Xn
T Kce ‰de Š :ˆ …Lg ‡ Qg ‰de Š† C Lg ‡ Qg ‰de Š ‡ Gg ‰ng Š ‡ STg DSg ˆ Kcg ‰de Š g

…57†

g

and, for the mixed one, h Xn
T Kme ‰Ue Š ˆ dde Gg ‰ng Šd_ e ‡ ddTe Lg ‡ Qg ‰de Š n_ g ‡ ddTe …STg DSg †d_ e ‡ dnTg …Lg ‡ Qg ‰de Š†d_ e

Fn_ g

io

e;g

…58† which provides Kme ‰Ue Š

:ˆ

X g

F …Lg ‡ Qg ‰de Š†T

Lg ‡ Qg ‰de Š Gg ‰ng Š ‡ STg DSg

 ˆ

X g

Kmg ‰Ug Š:

…59†

The simple, constant piecewise interpolation that is used for the in plane stress ®eld gives us the possibility of performing the elimination of the stress variables (by static condensation) directly at the single Gauss point level.

3384

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

By static condensation of the columns and rows of matrix Kmg ‰Ug Š associated to stresses ng and by using the identity C ˆ F 1 , we derive the reduced, pseudo-compatible matrix which formally coincides with matrix Kcg ‰de Š de®ned in Eq. (57) in the compatible formulation (using the same numerical integration scheme with ng ˆ Cg ). 4.4.3. The secondary forces vector ^ e; U _ e Š (i.e the vector form of the third The element contribution to the secondary internal forces se ‰U variation gived by Eq. (11)) is de®ned by the energy equivalence X X ^ e; U _ e Š :ˆ ^ g; U _ g Š: _ :ˆ U000 ‰uŠ^ uudu dUTe se ‰U dUTg sg ‰U e

e;g

By some manipulation, we obtain for the compatible model n o X X
T sce ‰^ de ; d_ e Š ˆ ddTe GTg ‰C^e Šd_ e ‡ Qg ‰^ de ŠT C_ e ‡ Lg ‡ Qg ‰de Š C^_ g Ag ˆ ddTe scg ‰^de ; d_ e Š

…60†

and for the mixed one  o Xn  X ^ e; U _ eŠ ˆ ^ g; U _ g Š; sme ‰U ddTe GTg ‰^ ng Šd_ e ‡ QTg ‰^ de Šn_ g ‡ dnTg ^_ g Ag ˆ dUTg smg ‰U

…61†

e;g

e;g

e;g

i.e.

( ^ g; U _ g Š :ˆ smg ‰U

e;g

^_ g T Gg ‰^ ng Šd_ e ‡ QTg ‰^ de Šn_ g

) Ag :

…62†

4.5. Element assembling The element assembling and all operations related to the macro-element management are the same as already discussed in [16,37] to which the reader is referred for details. Note, however, that the rotation matrix Re introduced in [16] to de®ne the displacement ®eld transformation from the local element coordinate system to the global one, now has a di€erent appearance due to the introduction of the stress ne into the representation of Ue   0 I RM :ˆ n : 0 RC In is here a 12  12 identity matrix that allows the components of ne to remain unchanged, and RC is the rotation matrix for compatible formulation already de®ned in [16,37]. 4.6. Pseudo-compatible solution scheme It is well known that a mixed problem can, conveniently, be solved by using a pseudo-compatible strategy in order to avoid the computational extra costs due to the increase in the number of variables [23,25]. This solution scheme allows the use of the compatible matrix, which is banded and positive de®nite for the solution of the linear systems and the eigenvalue problem generated by the perturbation algorithm.We now summarize this strategy in the solution of a system of equations: K‰kŠu ˆ f;

f :ˆ ff n ; f d g;

u :ˆ fn; dg

using a Newton±Raphson iteration scheme. We obtain: K0 u_ j ˆ rj ;

rj :ˆ f

sj ;

uj‡1 ˆ uj ‡ u_ j ;

…63†

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

3385

where j is the iteration index, K0 :ˆ K‰k0 Š and #   X" Fng;j ‡ …Lg ‡ Qg;j †dj sj;n sj :ˆ K‰kŠuj ˆ ˆ T sj;d …Lg ‡ Qg;j † ng;j ‡ …G‰ng;j Š ‡ STg DSg †dj;e e;g is the incremental force vector associated to uj (sj;n associated to stress variables and sj;d to displacement ones). Without any loss of generality we can refer to the case of a single panel with coincident local and global coordinates system (that is R ˆ I). The ®rst of Eq. (63) is     X F Lg ‡ Qg;0 n_ g;j f n sj;n ˆ : …64† T f d sj;d …Lg ‡ Qg;0 † …Gg ‰ng Š ‡ STg DSg † d_ e;j e;g After elimination of the stress variables, we obtain X e;g

F T …Lg ‡ Qg;0 †

Lg ‡ Qg;0 Kcg;0



n_ g;j d_ e;j

 ˆ

X e;g

f g;d

f g;n ‡ Fng;j …Lg ‡ Qg;j †dj T ‡ …Lg ‡ Qg;0 † Cf g;n ‡ …Qg;0 Qg;j †ng;j

 sg;j

;

where sg;j :ˆ …Lg ‡ Qg;0 †T C…Lg ‡ Qg;j †de;j ‡ …G‰ng;j Š ‡ STg DSg †dj;e :

…65†

After performing the static condensation at the element-level we get the pseudo-compatible matrix described only in terms of the kinematical ones. Then at the global-level we can proceed as usual by assembling the element contributions relative to the displacement degrees of freedom (i.e. the sti€ness matrix structure is the same as the compatible one). Again, as in the compatible approach, the global solution is expressed by using only the kinematical variables. It is possible then to proceed to the evaluation of the inplane stress ®eld in each Gauss point, by back-substitution of the displacement solution at the Gauss point level. i.e. h i ng;j‡1 ˆ ng;j ‡ n_ g;j ˆ C Lg …dj ‡ d_ j † ‡ Qg;0 d_ j ‡ Qg ‰dj Šdj …66† Cf g;n : For comparison the compatible model stresses solution is given by:     1 1 ncj‡1 ˆ Cj‡1 ˆ C Lg ‡ Qg ‰dj‡1 Š dj‡1 ˆ C Lg ‡ Qg ‰dj ‡ d_ j Š …dj ‡ d_ j †; 2 2 _ dj‡1 ˆ dj ‡ dj : The convenience of the pseudo-compatible implementation of the mixed formulation is evident. In fact with only small changes it is possible to convert a code which is based on the compatible approach to the described mixed extrapolation approach. 4.7. Compatible formulation with frozen geometry We also brie¯y recall the formal expressions belonging to the frozen geometry perturbation analysis, that are used to cure the extrapolation locking e€ect in compatible FE formulation [16,25]. Due to the assumption of the neglect of prebuckling displacement the quantity (47) becomes _ g ˆ Lg d_ e and the sti€ness matrix assumes the form o Xn Kfc LTg CLg ‡ Gg ‰ng Š ‡ STg DSg e ‰de Š :ˆ g

3386

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

while for the structural response vector we have o Xn s_ fc STg Dv_ g ‡ Gg ‰Cg Šd_ e ‡ LTg CLg d_ e : e ‰de Š :ˆ g

Further simpli®cations can be obtained in the evaluation of the third variation along the fundamental path [16]. The axial stresses being evaluated as: ncj‡1 ˆ Cj‡1 ˆ CLg …dj ‡ d_ j †:

…67†

5. Numerical results In this section we show the accuracy of the mixed formulation in comparison to the frozen con®guration method by performing some numerical tests. We start with simple tests simulating plane frame structures. In such a simpli®ed context it is easy to obtain the `exact' reference solutions either by the arc-length method [23] or by the mixed formulation of the perturbation approach for a kinematical exact two dimensional beam model [25,33]. In such a simpli®ed framework, the accuracy of the mixed approach will be veri®ed for the analysis of structures with a nonlinear precritical behavior. After these examples, more complex panel structures will be analyzed and the results compared with solutions obtained by commercial codes such as MSC±NASTRAN. In order to simplify the notation, we denote with Frozen Con®guration (FC) and Mixed (M) the results obtained with the perturbative Koiter Analysis of Slender Panels (KASP code) and with NASTRAN the results obtained using MSC±NASTRAN. The results obtained with the beam model implementation of mixed and frozen con®guration perturbative scheme [25] are labeled with FC±Beam and M±Beam, respectively. We label with Riks the path-following beam model results [23]. For comparison, we also report the results of the buckling loads obtained with NASTRAN. It is worth noting that the NASTRAN buckling loads are obtained with a linearized buckling analysis, i.e. with the implicit assumption of neglecting prebuckling displacements. 5.1. Euler beam with membranal constraints The geometry, constraints and load conditions for the ®rst numerical test are presented in Fig. 7. The structure, is a classical Euler beam with built-in edges for which the presence of the inclined lateral beams lead to the desired condition of stresses redundancy. Thus, we are able to evaluate the proposed KASP implementation in a context in which the error due to kinematical incoherences has no in¯uence on the solution. It is then possible to compare the results using an `exact' two-dimensional beam model solution both for perturbation and for path-following analysis.

Fig. 7. A simple frame with the ®rst buckling mode.

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

3387

The rectangular cross-section of the beams has a width b ˆ 10:0 while the thickness is t ˆ 2:0 for the inclined beams and t ˆ 1:0 for the central one. The structure consists of an isotropic homogeneous material with Young modulus E ˆ 1000. To simulate the two-dimensional behavior a Poisson ratio m ˆ 0:0 is assumed. Three di€erent meshes are considered: the ®rst grid labeled with G1 is given from a partition of …5 ‡ 10 ‡ 5†  2 elements (10 elements are used in the central beam), the second grid labeled G2 is built of …5 ‡ 20 ‡ 5†  2 partition and the third one, labeled G3, is given by …5 ‡ 40 ‡ 5†  2 elements. The solution is not dependent on the transverse direction therefore only two elements are needed in this direction (that is the minimum number required by the current KASP implementation of the HC interpolation). The analogous two-dimensional structure is built up of 2 ‡ 10 ‡ 2 beam elements. In Fig. 7 the ®rst buckling mode is plotted while in Table 1 some of the energy quantities relevant for the perturbation method are listed. Due to the linear precritical behavior of the structure the di€erences from the beam solution are negligible for both the mixed and the frozen geometry solution strategy. All meshes show an irrelevant error with respect to the reference, plane beam, solution. The perfect agreement of the equilibrium path obtained with both the FC and M solution schemes with the `exact' reference solution given by the Riks code is highlighted in Fig. 8 (a horizontal load imperfection with k~ p ˆ k0:001 is assumed). 5.2. Test2: a shallow arc The shallow arc depicted in Fig. 9 is a good test to illustrate the accuracy of the proposed perturbation method for structures characterized by a strong nonlinear behavior. The same test was analyzed in [25] by means of a plane beam element. The adopted geometrical ratio h=l (equal to 0:015) makes the pre-critical behavior of the structures strongly nonlinear. Table 1 Euler beam: comparison of relevant perturbative energy quantities

kb A001 A111 B0011 B0111 B1111

M-beam

FC-beam

M-G1

M-G2

M-G3

FC-G3

7.9433 0.0 0.0 4.410 0.0 1024.65

7.9433 0.0 0.0 4.410 0.0 1024.79

8.2080 0 0 4.510 0 1024.95

8.0089 0 0 4.510 0 1024.35

7.9597 0 0 4.510 0 1024.34

7.9597 0 0 4.510 0 1024.51

5

5

5

5

5

Fig. 8. Equilibrium path for Euler beam with membranal constraints (k~ p ˆ 0:001k).

5

3388

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

Fig. 9. Shallow arc with buckling modes.

Table 2 Shallow arc: comparison of relevant energy perturbation quantities

k1 k2 A001 A111 B0011 B0111 B1111

M-beam

FC-beam

M-G1

M-G2

M-G3

FC-G3

22.0202 30.1738 0.01934 48.4907 )6.26 10 3.51 10 8.52431

30.1738 46.6956 0.01838 177.138 0 0 58.0400

22.0763 30.8943 0.01936 48.5102 6.16 10 0 8.15312

21.9962 30.3504 0.01935 48.4819 6.30 10 0 8.43718

21.9752 30.2122 0.01934 48.4745 6.34 10 0 8.51369

30.2122 46.7597 0.01841 177.655 0 0 60.7125

4 4

4

4

4

As a consequence the results provided by the frozen and mixed formulations are completely di€erent. In Table 2 a more detailed comparison is reported showing the energetic quantities evaluated by the di€erent methods and by using di€erent meshes. In particular the G1 mesh is given by …6 ‡ 6†  2 HC ®nite elements, the G2 mesh by …12 ‡ 12†  2 HC ®nite elements, the mesh labeled G3 corresponds to a …25 ‡ 25†  2 grid. A mesh of 6 ‡ 6 beam elements was adopted for the analysis performed by the twodimensional model. Note that the frozen geometry model detects two buckling modes: the ®rst (skew-symmetric) at kb  30 and the second (symmetric) at kb  46. In contrast mixed formulation obtains kb  22 for the ®rst (symmetric) critical load. The frozen approach overestimates the limit-load value klim , as reported in Fig. 10. The same ®gure shows the good agreement between the results obtained by the KASP mixed model and those obtained with the beam model with Riks or perturbation analysis. In particular for the limit-load we obtain the same value klim  11:20. The results related to the FC approach, listed in Table 2, are evaluated for the symmetric (second) mode. 5.3. Angle section beam The simple structures analyzed in Section 5.2 were chosen in order to show the advantages of the proposed mixed method when standard reference solutions are available. Now test cases concerning applications of interest in the present work are considered. The ®rst test-case is an angle section beam whose geometry and load conditions are shown in Fig. 11 and are similar to the structure analysed in [38]. To magnify the e€ect of the nonlinear precritical behavior, due to the applied load eccentricity, the panels are 2:5 thick with respect to the value of 6:22 used in [38].

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

3389

Fig. 10. Equilibrium paths. Test 2.

Fig. 11. Angle section beam: geometry and buckling mode.

According to the latter thickness value the structure presents an unstable postcritical behavior that makes the di€erence (in the evaluation of the limit-load) between frozen and mixed approaches more evident. Because of the particular oblong form, required to highlight the e€ect of the precritical displacements and because of the presence of a global buckling mode, the structure could not be studied with the Simpli®ed-Lagrangian plate model due to the absence of the stress redistribution effect. This problem has been overcome by using a technical Complete-Lagrangian strain measure which is exact when the overall structural behavior is characterized essentially by displacements in the plane of the panels [16]. The results, see Fig. 12, con®rm the good behavior of the M strategy with respect to the FC one if compared with the equilibrium path obtained by NASTRAN code, which is assumed to be the exact reference. In Table 3 the most relevant quantities used in the perturbation solution are listed both for M and FC strategies using two di€erent meshes: the G1 is a …3 ‡ 5†  30 HC ®nite element grid, the G2 is a …4 ‡ 6†  60 HC ®nite element grid. The NASTRAN code results are obtained with a grid of …4 ‡ 6†  60 CQUAD4 elements. The mixed formulation is clearly more accurate than the frozen formulation.

3390

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

Fig. 12. Angle section beam: equilibrium paths obtained by the mixed and frozen model and by NASTRAN code.

Table 3 Angle section beam: comparison of relevant perturbative energy quantities M-G1

M-G2

FC-G1

FC-G2

NASTRAN

8722.7 8063.3 1.5258 19.388 0.0 15.908

8720.0 8053.6 1.5227 19.345 0.0 16.354

8716.9 8280.6 1.5274 16.426 0.0 5.5936

8714.1 8193.7 1.5287 16.418 0.0 8.2509

8757.4 8038.3 ± ± ± ±

kb klim A001  104 A111 B0011 B1111

5.4. Channel section The last test is a channel section, whose geometrical and mechanical characteristics are depicted in Fig. 13, subjected to a k-increasing compressed axial load distribution such that p^ ˆ 125. The structure is analyzed with three di€erent load imperfections (see Fig. 13) in order to investigate the performances of the various strategies when the initial precritical behavior addresses the direction of the torsional, local or ¯exural modes, respectively. The imperfection loads for the three cases analyzed consist of a system of two forces (see Fig. 13) applied to points A and B of the section at x ˆ l=2  550. For the torsional imperfection

Fig. 13. Channel section and force imperfection.

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

3391

we make p~A ˆ p~B ˆ k acting along the z direction, for the local one p~A ˆ p~B ˆ 2k acting along the y direction and for the ¯exural imperfection p~A ˆ p~B ˆ 2k acting along the z direction. In order to obtain a nonlinear precritical behavior and then highlight the better performance of the proposed method, the applied imperfection loads were considered together with the axial ones instead of in the postprocessing phase as would be convenient in a perturbation solution. As in the previous test the results obtained with M and FC strategies are compared with the NASTRAN solution. For the analysis performed with the KASP code a mesh of …3 ‡ 9 ‡ 3†  55 HC ®nite elements was used, for the NASTRAN analysis a …3 ‡ 9 ‡ 3†  54 grid of CQUAD4 elements was adopted. 5.4.1. Case A: torsional imperfection The ®rst four buckling modes obtained by the mixed formulation are drafted in Fig. 14. The buckling loads and a description of the the mode shape obtained with the three di€erent codes are listed in Table 4. The values and shapes of the buckling modes obtained with NASTRAN, M and FC are similar and the differences are essentially due to the different ®nite element implementations. It is worth noting that NASTRAN code uses a linearized buckling analysis neglecting the prebuckling deformations. In this case the results of the buckling analysis must be compared with those of the FC ones. The differences between M and FC are due to the effect of the nonlinear precritical behavior. The accuracy of the mixed strategy in the evaluation of the limit-load and of the initial postcritical behavior is highlighted in Figs. 15 and 16. In particular the evaluation of the limit-load of the proposed formulation agrees perfectly with the NASTRAN values (see Table 5) proving the reliability of the mixed approach. The equilibrium paths, obtained by both M and FC strategies are based on an expansion of the ®rst two or the ®rst four buckling modes. As highlighted in Figs. 15 and 16 the perturbation approaches show the strong e€ect of modal interaction between the third local mode and the ®rst two ¯exural±torsional

Fig. 14. Channel section: buckling modes for the torsional imperfection.

Table 4 Buckling values for the channel section with torsional imperfection Mode

v_ 1 v_ 2 v_ 3 v_ 4

Mode type

kcrit: M

FC

NASTRAN

1264.9 1833.9 2808.5 2823.9

1289.9 1711.8 2912.1 2930.0

1296.0 1715.6 2842.4 2863.5

Flexural Torsional Local (7 half-waves) Local (8 half-waves)

3392

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

Fig. 15. Case A: equilibrium path k

wA , k

wB .

Fig. 16. Case A: equilibrium path k

vB , k

uC .

Table 5 Limit-load values for the channel section with torsional imperfection klim M

FC

NASTRAN

1128:67

1204.17

1116.96

global modes (see Fig. 17). This mode interaction leads to unstable behavior. The NASTRAN equilibrium path has a di€erent evolution in the postbuckling range in uC and w while the k vB equilibrium path displays a good agreement also in the post-buckling range. These differences could be due to the presence of high order terms neglected in the perturbation analysis.

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

Fig. 17. Case A: equilibrium path (n1

n2 ), (n1

3393

n3 ).

5.4.2. Case B: local imperfection In Fig. 18 the ®rst eight buckling modes of the mixed formulation are shown. In this case there are also other modes for slightly higher k values. Table 6 reports the buckling load values and a description of the mode shape obtained with the three di€erent codes. There are minimal di€erences in the buckling loads obtained with the three approaches. Comparing the eigenvalues of NASTRAN with those of the FC ones there is a good agreement while the di€erences, related to the local modes, are due to the dissimilar ®nite element implementation. Small di€erences in the eigenvalues obtained by M and FC, due to the e€ect of the nonlinear precritical behavior, can be observed.

Fig. 18. Channel section: buckling modes for the local imperfection.

3394

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

Table 6 Buckling values for the channel section with local imperfection Mode

v_ 1 v_ 2 v_ 3 v_ 4 v_ 5 v_ 6 v_ 7 v_ 8

Mode type

kcrit: M

FC

NASTRAN

1293.2 1409.5 3150.7 3150.9 3206.7 3213.5 3278.4 3286.1

1289.9 1387.7 3134.3 3148.6 3161.1 3175.9 3213.8 3242.9

1296.1 1391.2 3084.8 3097.4 3112.9 3124.3 3167.9 3187.1

Flexural Torsional Local (13 half-waves) Local (14 half-waves) Local (12 half-waves) Local (11 half-waves) Local (15 half-waves) Local (10 half-waves)

As highlighted by the equilibrium paths (see Figs. 19 and 20) obtained taking into account the ®rst two or all the eight modes there is a strong e€ect of modal interaction between the ®rst ¯exural global mode and the third, sixth and seventh local modes, as is shown in Fig. 20. This interaction makes the equilibrium path unstable. There is a good agreement between the M and the NASTRAN results in the evaluation of the limit-load (see also Table 7) while the FC results are clearly rougher. The evaluation of the initial postcritical behavior obtained with the perturbation analysis detects interaction between local and global modes at a point of the equilibrium path beyond the point detected by NASTRAN. This e€ect can be attributed to the restricted subspace used with the perturbation approach, as can be observed when the subspace is increased to 16 buckling modes. It is worth noting that, when the number of buckling modes increase, the time required for the evaluation of the third- and fourth-order terms Aijk and Bijhk is a signi®cant percentage of the overall time analysis. In the context of simpli®ed kinematic relationships the mixed formulation gives an additional advantage due to the absence of the fourth variation terms in Bijhk . Although an exact comparison of the computational time obtained with the M and FC schemes is not possible due to some di€erences in code implementation, these computational time advantages are highlighted in Table 8 in which the single step time analysis is presented both for M and FC strategies.

Fig. 19. Case B: equilibrium path k

wA .

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

Fig. 20. Case B: equilibrium path k

uC , (n1

n3 ), (n1

n6 ), (n1

3395

n7 ).

Table 7 Limit-load values for channel section with local imperfection klim M

FC

NASTRAN

1246:1

1284.8

1235.42

Table 8 Time analysis for channel section with local imperfection CPU-time 8-mode

u^ v_ i wij Aijk Bijhk

16-mode

M

FC

M

FC

5 175 59 1 6

4 230 41 2 10

5 323 358 6 75

4 308 242 18 121

It is important to note how the di€erences in the computational time required by the eigenvalues analysis are due to the di€erent representation of the critical subspace in which the subspace iteration is performed. The greater computational time required for the evaluation of the orthogonal displacement wij is due both to the evaluation of the term w00 , not evaluated by the FC code, and to the heavier evaluation of the energy expressions. The di€erences in computational time required for the evaluation of the cubic terms is related to some implementational di€erences (FC code needs to recognize which expression has to be used for the third variation [16]). Finally in the evaluation of the fourth-order terms, the mixed approach has to evaluate all terms B00ij , identical to zero in the frozen one. Notwithstanding this the computational time is almost 60% lower than with the FC code.

3396

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

Fig. 21. Channel section: buckling modes for the ¯exural imperfection.

Table 9 Buckling values for channel section with ¯exural imperfection Mode

v_ 1 v_ 2 v_ 3 v_ 4

kcrit:

Mode type

Mixed

Frozen

NASTRAN

1291.6 1396.0 1994.4 2046.4

1292.2 1192.5 1992.8 2044.7

1298.2 1195.2 1919.4 2094.5

Flexural Torsional Local (3 half-waves) Local (4 half-waves)

5.4.3. Case C: ¯exural imperfection In Fig. 21 the ®rst four buckling modes detected by mixed formulation are depicted. In Table 9 the buckling loads and a description of the mode shape obtained with the three di€erent codes are reported. Note the agreement between FC and NASTRAN results with respect to the buckling loads, the differences on the local modes being essentially due to the di€erent ®nite element technologies. A greater di€erence with respect to the M approach, that detects as the ®rst mode the torsional one (the second for FC and MSC±NASTRAN (see Table 9)) can also be observed. There is a good agreement between M and NASTRAN in the evaluation of the limit-load as shown in Table 10. The equilibrium paths (see Figs. 22 and 23) show the best approximation given by the mixed formulation with respect to the frozen one. In this case, however, the di€erences are of minor relevance. The e€ect of neglecting prebuckling deformation is highlighted in Fig. 23 where the solutions near the bifurcation point are ampli®ed. The postcritical behavior is characterized by the interaction between the ®rst (¯exural) and the third (local) critical modes (see Fig. 24). The results show a good agreement with those obtained with NASTRAN especially in the vB displacement component. A few di€erences can be observed in the uC and wA components with the usual deeper descending path depicted by the perturbation analysis. Table 10 Limit-load values channel section with ¯exural imperfection kmax Mixed

NASTRAN

Frozen conf.

900.3

896.3

908.4

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

Fig. 22. Case C: equilibrium path k

Fig. 23. Case C: equilibrium path k

wA , k

vB .

uC , zoom on the prebuckling path.

Fig. 24. Case C: equilibrium path (n1

n3 ).

3397

3398

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

6. Conclusion The paper extends the mixed version of the perturbation algorithm presented in [25] to the structures already studied in [15,16], allowing a cure of locking extrapolation problems, typical of the compatible formulation, without using an ad hoc hypothesis such as the frozen con®guration assumption. The latter, already adopted in previous works [15,16] on the same structural type, even if useful for eliminating locking problems can give some unreliable results for structures characterized by nonlinear precritical behavior. The possibility of solving the mixed FE algebraic system and eigenvalues search, using a pseudo-compatible solution scheme, makes the analysis fast with no computational extra cost with respect to the compatible one. When used, on the basis of simpli®ed quadratic strain±displacement relationships, the mixed formulation has the further advantage of having all energy variations higher than third-order zeroed. This produces noticeable computational simpli®cation such as a natural linearization of the eigenvalue problem that is used for the buckling analysis. All the numerical tests performed show the reliability and accuracy of the mixed formulation with respect to the frozen one. While the limit-load and the initial postcritical behavior obtained with the proposed formulation always completely agrees with the path-following NASTRAN results, di€erences in the postcritical equilibrium path evolution can be observed with the perturbative results that exhibit a more unstable postbuckling behavior. To understand this deviation some investigations will have to be carried out. A ®rst step in this direction is the elimination of the FE implementation di€erences of the path-following versus perturbation codes (a path-following version of the KASP code is being developed). A rational kinematical plate model, suitable for use in a perturbation approach, is another important topic which is being investigated. Acknowledgements This research was supported by MURST and Brite±Euram APRICOS (EC BE95-1017). The author is very grateful to Ing. E. Feraco for his collaboration in developing the MSC±NASTRAN numerical tests. A special acknowledgement is also due to Prof. R. Casciaro, of the University of Calabria (Italy), for his precious observations and suggestions. References [1] W.T. Koiter, On the stability of elastic equilibrium, Thesis, Delft, 1945. English translation NASA TT-F10, 883 (1967) and AFFDL/TR70-25, 1970. [2] W.T. Koiter, Current trends in the theory of buckling, in: B. Budiansky (Ed.), Buckling of Structures, Proceedings of IUTAM Symposium, Cambridge, 17±21 June 1974, Springer, Berlin, 1976. [3] B. Budiansky, Theory of Buckling and Post-buckling of Elastic Structures, in: Advances in Applied Mechanics, vol. 14, Academic Press, New York, 1974. [4] E. Riks, An incremental approach to the solution of snapping and buckling problems, Int. J. Solids Struct. 15 (1979) 529±551. [5] E. Riks, Buckling analysis of elastic structures; a computational approach, in: E. van der Giessen, T.Y Wu (Eds.), Advances in Applied Mechanics, Academic Press, New York, 1997. [6] D. Ho, Higher order approximations in the calculation of elastic buckling loads of imperfect systems, Int. J. Non-Linear Mech. 6 (1971) 649±661. [7] D. Ho, The in¯uence of imperfections on systems with coincident buckling loads, Int. J. Non-linear Mech. 7 (1972) 311±321. [8] D. Ho, Buckling load of non-linear systems with multiple eigenvalues, Int. J. Solids Struct. 10 (1974) 1315±1330. [9] W.T. Koiter, Some properties of (completely) symmetric multilinear forms with an application to Ho's theorem for multi-mode buckling, lab. rep. no 587, Technologic University of Delft, The Netherlands, April 1976. [10] G. Salerno, R. Casciaro, Mode jumping and attrative paths in multimode elastic buckling, Int. J. Numer. Methods Engrg. 40 (1997) 833±861. [11] G. Salerno, Secondary bifurcations and attractive paths in monomodal reduction of coupled instability, in: Proceedings of the Second International Conference on Coupled Instabilities in Metal Structures, Liege, Belgium, September 1996, pp. 5±7. [12] G. Salerno, G.Colletta, R.Casciaro, Attractive post-critical paths and stochastic imperfection sensitivity analysis of nonlinear elastic structures within an FEM context, in: Proceedings of the Second ECCOMAS International Conference on Numerical Methods in Engineering, Paris, France, September 1996.

G. Garcea / Comput. Methods Appl. Mech. Engrg. 190 (2001) 3369±3399

3399

[13] R. Casciaro, A.D. Lanzo, G. Salerno, Computational problems in elastic structural stability, in: C. Carmignani, G. Maino (Eds.), Nonlinear Problems in Engineering, Proceedings ENEA Workshops on Nonlinear Dynamics, vol. 4, Roma, 6±7 May 1991, World Scienti®c, Singapore, 1991. [14] R. Casciaro, G. Salerno, A.D. Lanzo, Finite element asymptotic analysis of slender elastic structures: a simple approach, Int. J. Numer. Methods Engrg. 35 (1992) 1397±1426. [15] A.D. Lanzo, G. Garcea, R. Casciaro, Asymptotic post-buckling analysis of rectangular plates by HC ®nite elements, Int. J. Numer. Methods Engrg. 38 (1995) 2325±2345. [16] A.D. Lanzo, G. Garcea, Koiter's analysis of thin-walled structures by s ®nite element approach, Int. J. Numer. Methods Engrg. 39 (1996) 3007±3031. [17] A.D. Lanzo, G. Garcea, Asymptotic post-buckling analysis of thin-walled structures, in: Proceedings of the International Conference on Mechanics of Solids and Materials Engineering, Singapore, June 1995, pp. 5±7. [18] A.D. Lanzo, G. Garcea, Coupled instabilities in thin-walled structures by a FEM implementation of Koiter's approach, in: Proceedings of the Second International Conference on Coupled Instabilities in Metal Structures, Liege, Belgium, September 1996, pp. 5±7. [19] G. Salerno, A.D. Lanzo, Sull'insensibilita post-critica all'esatezza geometrica ed al modello FEM di strutture staticamente indeterminate, Congresso AIMETA '97, Siena, Sett. 1997. [20] M. Pignataro, A. Di Carlo, R. Casciaro, On nonlinear beam model from the point of view of computational post-buckling analysis, Int. J. Solids Struct. 18 (4) (1982) 327±347. [21] R. Casciaro, G. Garcea, G. Attanasio, F. Giordano, Perturbation approach to elastic post-buckling analysis, in: Proceedings of MSC European User Group, Rome, 1995, Comput. Struct. 66 (1998) 585±595. [22] F. Brezzi, M. Cornalba, A. Di Carlo, How to get around a simple quadratic fold, Numer. Math. 48 (1986) 417±427. [23] G. Garcea, G.A. Trun®o, R. Casciaro, Mixed formulation and locking in path following nonlinear analysis, Comput. Methods Appl. Mech. Engrg. 165 (1±4) (1998) 247±272. [24] G. Garcea, G. Salerno, Sanitizing of locking in Koiter perturbation analysis through mixed formulation, Fourth World Congress on Computational Mechanics, Buenos Aires, Argentina, 1998. [25] G. Garcea, G. Salerno, R. Casciaro, Extrapolation locking and its sanitization in Koiter's asymptotic analysis, Comput. Methods. Appl. Mech. Engrg. 180 (1-2) (1999) 137±167. [26] S. Sridharan, Doubly symmetric interactive buckling of plate structures, Int. J. Solids Struct. 19 (1983) 625±641. [27] R. Benito, S. Sridharan, Interactive buckling analysis with ®nite strip, Int. J. Numer. Methods Engrg. 21 (1985) 145±171. [28] M. Ali, S. Sridharan, A versatile model for interactive buckling of columns and beam-columns, Int. J. Solids Struct. 24 (1988) 481± 496. [29] P. Goltermann, H. Mùllmann, Interactive buckling in thin-walled beams ± I. Theory, II Application, Int. J. Solids Struct. 25 (1989) 715±749. [30] G.M. van Erp, Advanced buckling analyses of beam with arbitrary cross-section, Ph.D. Thesis, Eindhoven University of Technology, 1989. [31] G.M. van Erp, C.M. Menken, Initial post-buckling analysis with the spline ®nite-strip method, Comput. Struct. 40 (5) (1991) 1193±1201. [32] R. Casciaro, A.D. Lanzo, G. Salerno, A new algorithm for nonlinear symmetric eigenvalue problems, Rep. 122, Dipartimento di Strutture, Universita della Calabria, April 1990. [33] G. Salerno, A.D. Lanzo, A nonlinear beam ®nite element for the post-buckling analysis of plane frames by Koiter's perturbation approach, Comput. Methods Appl. Mech. Engrg. 146 (1997) 325±349. [34] E. Byskov, Smooth postbuckling stresses by a modi®ed ®nite element method, Int. J. Numer. Methods Engrg. 28 (1989) 2877± 2888. [35] A.K. Noor, Recent advances and application of reduction methods, Appl. Mech. Rev. 5, May 1994, 125±146. [36] M. Aristodemo, A high-continuity ®nite element model for two-dimensional elastic problems, Comput. Struct. 21 (1985) 987±993. [37] A. Bilotta, G. Garcea, G.A. Trun®o, R. Casciaro, Postcritical analysis of thin walled structures by KASP code, APRICOS±DE± 1.3±03±1/UNICAL, Universita della Calabria, Report 185, September 1997. [38] S.L. Chan, S. Kitipornchai, Geometric nonlinear analysis of asymmetric thin-walled beam-columns, Engrg. Struct. 9 (4) (1987) 243±254.