A perturbation approach for geometrically nonlinear structural analysis using a general purpose finite element code

(1)

A Perturbation Approach for Geometrically

Nonlinear Structural Analysis Using a

General Purpose Finite Element Code

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magniﬁcus Prof.dr.ir. J.T. Fokkema, voorzitter van het College voor Promoties,

in het openbaar te verdedigen op 21 december om 10.00 uur

door

Tanvir RAHMAN

Master of Science in Computational Mechanics of Materials and Structures, University of Stuttgart, Stuttgart, Germany

(2)

Prof.dr. Z. G¨urdal

Samenstelling promotiecommissie:

Rector Magniﬁcus, voorzitter

Prof.dr. Z. G¨urdal, Technische Universiteit Delft, promotor Dr.ir. E.L. Jansen, Technische Universiteit Delft, copromotor Prof.dr. J. Arbocz, Technische Universiteit Delft

Prof.dr.ir. A. van Keulen, Technische Universiteit Delft Prof.dr.ir. A. de Boer, Universiteit Twente

Prof.dr. H. Nijmeijer, Technische Universiteit Eindhoven Dipl.ing. A. Newerla, ESA-ESTEC, adviseur

Publisher: TUDelft, Faculteit Luchtvaart-en Ruimtevaarttechiek Printed by: Ipskamp Drukkers B.V., Enschede, The Netherlands ISBN 978-90-9024951-3

Keywords: Perturbation, Postbuckling, Geometric nonlinearity, Finite Elements, Thin-walled Structures.

This research is supported by the Dutch Technology Foundation STW, applied sci-ence division of NWO and the Technology Program of the Ministry of Economic Aﬀairs (STW project 06613).

Copyright c°2009 by Tanvir Rahman

All rights reserved. No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic or mechanical, in-cluding photocopying, recording or by any information storage and retrieval system, without written permission from the publisher: TUDelft, Faculteit Luchtvaart-en Ruimtevaarttechiek

(3)

(4)

(5)

Acknowledgment

At the completion of this thesis, ﬁrst of all, I would like to thank my supervi-sor Prof.dr. Z.G¨urdal for giving me the opportunity to conduct this research in his group and for his continuous guidance. Special thanks to my direct supervisor Dr. E.L. Jansen who guided me all through with his valuable advices and ideas, which often extended from the research topic to the family matters. My sincere gratitude to the committee members for serving in the examination committee and for reading and providing their useful comments on the thesis, Prof.dr. J. Arbocz, Prof.dr.ir. A. van Keulen, Prof.dr.ir. A. de Boer, Prof.dr. H. Nijmeijer, and espe-cially to my adviser in the committee, Dipl.ing. A. Newerla for his very detailed recommendations.

I am indebted to the Dutch Technology Foundation, STW for the ﬁnancial sup-port in this research and also to all the STW users committee members, Dr.ir. A. Bout, Ir. V. Bouwman, Dr. A.F. Calvi, Ir. J.P. Delsemme, Dr.ir. M.H. van Houten, Prof.dr. H. Nijmeijer, Dr.ir. G.J.M.A. Schreppers, Ir. H. Vasmel, Ir. J.J. Wijker, and to the program oﬃcer, Dr.ir. C.H.J. Meuleman for their valuable recommenda-tions through out the research project.

I would like to thank my former colleague Dr.ir. P. Tiso for his excellent earlier research work in our group which paved my way to a large extent for the present research. I would also like to thank Dr.ir. M.M. Abdalla and Ir. S.T. IJsselmuiden for their contributions in a part of this research. I would also like to extend my thanks to Ir. Wijtze Pieter Kikstra for explaining me the relevant parts of DIANA code whenever necessary.

Finally, I would like to thank Ms. A. van Lienden, Ms. A. de Gier, and Ms. L. Chant for their spontaneous help about administrative issues and all my colleagues, family members, relatives, and friends who always encouraged and inspired me.

(6)

(7)

Summary

In this thesis, a finite element based perturbation approach is presented for geomet-rically nonlinear analysis of thin-walled structures. Geometgeomet-rically nonlinear static and dynamic analyses are essential for this class of structures. Nowadays nonlin-ear analysis of thin-walled shell structures is often done using finite element based incremental-iterative procedures. However standard finite element based nonlinear analysis of many practical structures is still computationally expensive and not suit-able for repeated runs necessary for a design and optimization process, and also the interpretation of the results can be difficult. Koiter introduced a perturbation ap-proach where a perturbation expansion for the load and the displacement field are made around a known state (e.g. a bifurcation buckling point) to capture the es-sential nonlinear behavior of the structure. The final outcome is a computationally inexpensive reduced order model providing a good insight into the structural behav-ior in terms of the postbuckling a and b coefficients and the first order (the buckling modes) and the second order modes. Therefore, application of Koiters perturbation approach within a finite element framework can support standard nonlinear analysis by providing an effective solution to the aforementioned problems. The main focus of the present research is on the extensions which are necessary for this application The aim of the present research is to show the feasibility of implementing the perturbation approach in a general purpose finite element code in order to handle practical shell structures. In accordance with this perspective, the implementation of the present work is done in the development environment of a general purpose finite element code (DIANA) so that practical problems can be addressed. A general framework for the finite element implementation of Koiters perturbation approach is already available in DIANA. The present implementation makes use of that frame-work and does the necessary extensions. Two key aspects in this research are the inclusion of the effect of prebuckling nonlinearity and consideration of composite shells. Inclusion of the effect of prebuckling nonlinearity makes it possible to handle an important class of problems consisting of axially loaded cylindrical and conical shells where the prebuckling nonlinearity has a large influence on the buckling and postbuckling behavior. The implementation of the perturbation approach is done for a class of Mindlin type curved beam and shell elements and a Kirchhoff type triangular flat shell element.

The present thesis contains a part on static problems and a part on dynamic problems. In the part on static initial post-buckling problems, ﬁrst the perturba-tion method is explained within the single mode context (considering one buckling mode), including prebuckling nonlinearity and its ﬁnite element implementation.

(8)

The special case of a linear prebuckling state is also discussed. Various reference problems consisting of both isotropic and composite shells are considered. The postbuckling coefficients are compared with available results in literature and semi-analytical results. In all cases a reasonable agreement in the buckling and second order modes and the postbuckling coefficients is obtained. In case of the benchmark conical shell problem (representative of the 1/2 interstage of the Vega launcher) initial postbuckling response (based on single mode analysis) obtained using the present approach compares reasonably well with the full model analysis. A com-parison of computational time with full model analysis shows the significant gain achieved in computational cost with the present approach. The locking problem in the determination of postbuckling b coefficient is an important issue in the im-plementation of the perturbation approach in a finite element context. It is found that the Mindlin type shell and beam elements that are used in this research (where both in-plane/axial and out-of-plane/lateral displacements are interpolated inde-pendently) do not show the locking problem with respect to convergence of the b coefficient.

Next, in order to study modal interaction effects, the necessary extensions needed for multi mode analysis (in which more than one buckling mode is considered) are discussed. Considering a large number of modes in the perturbation analysis is computationally expensive. Therefore, a guideline is proposed for the selection of appropriate modes based on the minimum direction (the direction that leads to the steepest descent or the smallest ascent of the total potential energy in the buckling space) and uncoupled postbuckling coefficients (postbuckling coefficients associated with a single specific mode) and demonstrated using appropriate examples. Also, industrial benchmark problems consisting of a curved and a flat stiffened CFRP panel are considered. A reasonable agreement in the results is obtained in both cases. The efficiency of the reduced order model is also demonstrated by comparing the computational time with full model analysis.

Finally, as a specific application, axially loaded, variable stiffness composite pan-els are considered and the potential of the perturbation approach in the optimization process for the improvement of postbuckling performance of such panels is demon-strated. The expressions for in-plane postbuckling stiffness for plates under axial compression is derived both for the single mode and the multi mode case. In a specific multi mode case where clustering of buckling modes appears, the minimum direction or the worst direction (in the buckling space) of the structure are deter-mined for the initial postbuckling regime and the resulting postbuckling stiffness is computed. It is found that the panels considered having two nearly coincident buckling modes show the lowest postbuckling stiffness if they buckle in a shape with two half waves. This can be explained in terms of the second order stress distri-bution over the panel. It is an instance showing the type of insight one can gain about the structural behavior using the perturbation approach. It is also found that a compromise between buckling load and prebuckling stiffness may yield a better performance with respect to ultimate load (the failure load, in this case assumed to correspond with the onset of material degradation) and limit load (the maxi-mum load that the structure can encounter during its operating life) compared to

(9)

ix only buckling dominated design. Up to limit load level good agreement is obtained between the full model and the reduced order model analysis.

In the part of the thesis about dynamic problems, the perturbation approach developed earlier for static postbuckling problems is extended for handling dynamic problems by taking inertial terms into account. Dynamic buckling and nonlinear free vibration problems are considered. A reasonable agreement in the prediction of the dynamic buckling load is achieved when compared with full model transient analysis. The gain in computational cost is shown for a representative case. In case of the nonlinear vibration problems, second order vibration modes, the dynamic b coeﬃcients and the resulting backbone curves compare reasonably well with refer-ence results.

(10)

(11)

Samenvatting

In dit proefschrift wordt een storingsrekeningsaanpak binnen de eindige elementen methode gepresenteerd voor geometrisch niet-lineaire analyse van dunwandige con-structies. Geometrisch niet-lineaire statische en dynamische analyses zijn essentieel voor deze klasse van constructies. Tegenwoordig wordt niet-lineaire analyse van dun-wandige schaalconstructies vaak uitgevoerd via incrementeel-iteratieve procedures binnen de eindige elementen methode. Standaard niet-lineaire eindige elementen analyse van vele praktische constructies kost nog steeds veel rekentijd en is niet geschikt voor de herhaalde berekeningen die nodig zijn binnen een ontwerp- en opti-malisatieproces, en ook de interpretatie van de resultaten kan moeilijk zijn. Koiter introduceerde een storingsrekeningsaanpak waarbij een storingsexpansie voor de be-lasting en het verplaatsingsveld wordt gemaakt rond een bekende toestand (b.v. een bifurcatie-knik punt) om de essentie van het niet-lineaire gedrag te beschrijven. Het uiteindelijke resultaat is een wat rekentijd betreft goedkoop gereduceerde orde model dat een goed inzicht verschaft in het gedrag van de constructie in de vorm van de naknikcoefficienten a en b en de eerste orde velden (knikvormen) en tweede orde velden. Daarom kan de toepassing van Koiter’s storingsaanpak binnen een eindige elementen raamwerk de standaard niet-lineaire analyse ondersteunen door een effec-tieve oplossing te verschaffen voor de eerder vermelde problemen. De uitbreidingen die nodig zijn voor deze toepassing vormen de kern van het onderhavige onderzoek. Het doel van het huidige onderzoek is het aantonen van de realiseerbaarheid van het implementeren van de storingsaanpak in een ’general purpose’ eindige ele-menten programma om zo praktische schaalconstructies te kunnen aanpakken. Met het oog hierop worden de implementaties van het huidige werk gedaan binnen de on-twikkelomgeving van een ’general purpose’ eindige elementen programma (DIANA) zodat praktische problemen behandeld kunnen worden. Een algemeen raamwerk voor de eindige elementen implementatie van Koiter’s storingsaanpak is al beschik-baar binnen DIANA. De huidige implementatie maakt gebruik van dat raamwerk en zorgt voor de noodzakelijke aanpassingen. Twee essentiele aspecten in dit onder-zoek zijn het meenemen van het effect van niet-lineariteit in de grondtoestand en het beschouwen van composietschalen. Meenemen van het effect van de niet-lineariteit van de grondtoestand maakt het mogelijk om de belangrijke klasse van problemen te behandelen van axiaal belaste cilindrische en conische schalen, waarin de niet-lineariteit van de grondtoestand een grote invloed heeft op het knik- en naknikge-drag. De implementatie van de storingsaanpak wordt gedaan voor een klasse van gekromde balk- en schaalelementen van het Mindlin-type en een driehoekig vlak schaalelement van het Kirchhoff-type.

(12)

Het onderhavige proefschrift bevat een gedeelte over statische problemen en een gedeelte over dynamische problemen. In het gedeelte over statische initiele naknikproblemen wordt eerst de ’single mode’ aanpak (waarin een enkele knikvorm wordt beschouwd) van de storingsmethode uitgelegd, inclusief het meenemen van niet-lineariteit in de grondtoestand en de corresponderende eindige elementen im-plementatie. Het speciale geval van een lineaire grondtoestand wordt bediscussieerd. Verschillende referentieproblemen met zowel isotrope schalen als composietschalen worden bekeken. De naknikcoefficienten worden vergeleken met beschikbare resul-taten in de literatuur en met semi-analytische resulresul-taten. In alle gevallen wordt een redelijke overeenkomst in knikvormen, tweede orde velden, en naknikcoefficienten verkregen. In het geval van een ’benchmark’ probleem van een conische schaal (rep-resentatief voor de 1/2 tussentrap van het Vega lanceervoertuig) komt de initiele naknikresponsie (gebaseerd op ’single mode’ analyse) verkregen met de huidige aan-pak redelijk goed overeen met de analyse met het volledige model. Een vergelijking van de rekenkosten met die van de analyse met het volledige model toont de belan-grijke winst in rekentijd die kan worden bereikt bij de huidige aanpak. Het ’locking’ probleem in de bepaling van de naknikcoefficient b is een belangrijk aandachtspunt bij de implementatie van de storingsaanpak binnen een eindige elementen omgeving. Het blijkt dat de schaal- en balkelementen van het Mindlin type die in dit onder-zoek gebruikt worden (waarbij zowel de in-het-vlak / axiale als de uit-het-vlak / laterale verplaatsingen onafhankelijk geinterpoleerd worden) het locking probleem met betrekking tot het convergeren van de b coefficient niet vertonen.

Vervolgens worden de voor het bestuderen van modale interactie effecten noodza-kelijke uitbreidingen bediscussieerd voor de ’multi mode’ analyse (waarin meerdere knikmodes in beschouwing worden genomen). Het meenemen van een groot aantal knikvormen in de storingsanalyse kost veel rekentijd. Daarom wordt een richtlijn voorgesteld voor de selectie van geschikte knikvormen gebaseerd op de ’minimum richting’ (de richting die leidt tot de steilste daling of kleinste steiging van de totale potentiele energie in de ruimte van de knikvormen) en ongekoppelde naknikcoeffi-cienten (naknikcoeffinaknikcoeffi-cienten corresponderend met een enkele specifieke knikvorm) en geillustreerd met toepasselijke voorbeelden. Ook worden industriele ’benchmark’ problemen bestaand uit een gekromd en een vlak verstijfd CFRP paneel bekeken. Een redelijke overeenkomst tussen de resultaten wordt in beide gevallen verkregen. De efficientie van het gereduceerde orde model wordt ook gedemonstreerd door het vergelijken van de rekentijd met die van de analyse met het volledige model.

Tenslotte worden als een speciﬁeke toepassing in-het-vlak belaste, variabele sti-jfheid composiet panelen bekeken en wordt de potentie van de storingsaanpak in het optimalisatieproces voor het verbeteren van de prestatie in het naknikgebied van zulke panelen aangetoond. De uitdrukkingen voor de in-het-vlak postbuckling stijfheid voor platen onder in-het-vlak compressie wordt afgeleid zowel voor het ’single mode’ geval als voor het ’multi mode’ geval. In een speciﬁek multi mode geval waarin clusteren van knikmodes optreedt, worden de ’minimum richting’ of de ’slechtste richting’ (in de ruimte van knikmodes) van de constructie bepaald voor het initiele naknikgebied en de resulterende naknikstijfheid wordt berekend. Het blijkt dat de beschouwde panelen bij twee bijna samenvallende kniklasten de laagste

(13)

xiii naknikstijfheid hebben wanneer zij knikken in een knikvorm met twee halve golven. Dit kan verklaard met behulp van een beschouwing van de spanningsverdeling van de tweede orde velden. Dit is een voorbeeld van het type inzicht dat kan worden verkregen over het gedrag van de constructie bij het gebruik van de storingsaan-pak. Verder wordt gevonden dat een compromis tussen de kniklast en de stijfheid in de grondtoestand een betere prestatie met betrekking tot de ’ultimate load’ (de bezwijkbelasting, in dit geval aangenomen te corresponderen met de aanvang van materiaaldegradatie) en de ’limit load’ (de maximale belasting die de constructie in bedrijf kan tegenkomen) zou kunnen opleveren. Tot ’limit load’ level wordt een goede overeenkomst tussen het volledige model en het gereduceerde orde model verkregen. In het deel van het proefschrift over dynamische problemen wordt de storingsaan-pak eerder ontwikkeld voor statische naknikproblemen, uitgebreid voor de behan-deling van dynamische problemen door traagheidstermen in de beschouwing te be-trekken. Dynamische knikproblemen en niet-lineaire vrije trillingsproblemen worden bekeken. Een redelijke overeenkomst in de voorspelling van de dynamische kniklast wordt verkregen bij een vergelijking met een transiente analyse met het volledige model. De winst in rekentijd wordt getoond voor een representatief geval. In het geval van het niet-lineaire trillingsprobleem komen de tweede orde trillingsvormen en de dynamische b coeﬃcienten goed overeen met referentieresultaten.

(14)

(15)

Nomenclature

u, v Generalized displacements

σ, ² Generalized stress and strain

m The number of buckling modes or vibration modes considered

i, j, k, l Indices running from 1 to m,

subject to Einstein’s summation convention

I, J, K, L Indices running from 1 to m,

not subject to Einstein’s summation convention unless otherwise stated

u0, ²0, σ0 Prebuckling displacement, strain, and stress ﬁelds

ui, ²i, σi Buckling or ﬁrst order displacement, strain, and stress ﬁelds

uij, ²ij, σij Second order displacement, strain, and stress ﬁelds

H Linear operator used in the stress-strain relation;

also used to denote the height of cylindrical and conical shells

M Mass operator

t Time; also used to denote shell thickness

L1 Linear operator used in the linear part

of the strain-displacement relation

L11, L2 Quadratic operators used in the nonlinear part

of the strain-displacement relation

ξi Perturbation parameter associated

with buckling or vibration mode i

(20)

f External load vector

λ Load parameter, such that f = λf0

λc Bifurcation buckling load of a perfect structure

λd Dynamic buckling load of a perfect or

an imperfect structure

λs Static buckling load of an imperfect structure

λb Base state load level, where λb / λc

( )c, ( )b Prebuckling quantities evaluated respectively,

at λ = λc and λ = λb

g Right hand side force vector for

determination of the second order modes

φ A parameter used for determination of the second order mode, where φ/ 1

˙

( ) = _∂λ∂ ( ) in case of static postbuckling analysis ˙

( ) = _∂t∂( ) in case of dynamic buckling and vibration analyses

aijk, bijkl Static postbuckling coeﬃcients

aDijk, bDijkl Dynamic a and b coeﬃcients

αi, βi Imperfection form factors

P Total potential energy

Aijk, Aijkl Coeﬃcients in the expansion of P

ˆ

u Geometric imperfection shape ¯

ξ Imperfection amplitude

ei Minimum direction in the space of buckling modes

(21)

xxi

ω Circular frequency of vibration (large amplitude); also used for small amplitude vibrations

when the context is not related to nonlinear vibration problems

τ Non-dimensional time, where τ = ωt ˆ

u1 Spatial part of the vibration mode,

such that u1 = ˆu1cos τ

ˆ

u21, ˆu22 Spatial part of the second order mode

in case of nonlinear vibration problems such that u2 = ˆu21 + ˆu22cos 2τ

q Displacement in discretized form (nodal displacement vector)

Kxx, Kyy, Kzz, Matrices containing derivatives of the

Kxy, Kyz, Kzx interpolation polynomial functions

(not to be confused with stiﬀness matrices)

BL Linear part the of strain-displacement matrix

BN L Nonlinear part the of strain-displacement matrix

H Stress-strain matrix, such that σ = H²

KMe, KM Element and assembled material stiffness matrices KDe, KD Element and assembled initial displacement matrices KGe, KG Element and assembled geometric stiffness matrices Kte, Kt Element and assembled tangent stiffness matrices KLde, KLd Element and assembled load stiffness matrices Me, M Element and assembled mass matrices

E11, E22 Young’s moduli of an orthotropic layer

G12 Shear modulus of an orthotropic layer

(22)

Xt Yield strength of an orthotropic material

in tension along X axis

Xc Yield strength of an orthotropic material

in compression along X axis

Yt Yield strength of an orthotropic material

in tension along Y axis

Yc Yield strength of an orthotropic material

in compression along Y axis

S Yield strength of an orthotropic material in shear

r Failure index

L Length of a cylindrical shell

R Radius of a cylindrical shell

c =p3(1− ν2

12)

θ Fiber orientation angle

Pcr Buckling pressure (force/unit area)

Ncl Bucking load (force/unit length)

N Number of circumferential waves in the buckling mode of cylindrical and conical shells

Nm Mesh size in terms of number of

divisions along a line or a curve

αs Semi-vertex angle of a conical shell

a, b Length and width of a rectangular plate, respectively

µ = a/b

h Thickness of a rectangular plate

ρ Density

(23)

xxiii

Kpr Prebuckling stiﬀness of an axially loaded panel

Kpo Postbuckling stiﬀness of an axially loaded panel

(24)

(25)

Chapter 1 Introduction

1.1 Background

Thin-walled structures such as plates and shells possess a high strength-over-weight and stiﬀness-over-weight ratio, and therefore are used as the primary components in weight critical structural applications in Aerospace, Mechanical, Marine, and Civil engineering. These structures are prone to static and dynamic buckling instabilities, and they can easily be excited into resonance. Often at the onset of buckling the stress level remains much lower than the yield stress and buckling is considered to be the key design criteria. Moreover, shell type structures exhibit unstable postbuckling behavior which makes them highly sensitive to small geometric or load imperfections. When loaded, thin walled structures may be subject to high out-of-plane displace-ments (compared to wall thickness) and cause geometrically nonlinear structural response. Therefore, geometrically nonlinear static and dynamic analyses are essen-tial for this class of structures. Important analysis cases in this ﬁeld include static buckling and postbuckling, dynamic buckling, nonlinear vibration, mode jumping (analyzed with a hybrid static-dynamic approach) and general transient analysis.

Nowadays nonlinear analysis of thin-walled shell structures are often done us-ing finite element based incremental-iterative procedures. However standard finite element based nonlinear analysis of many practical structures is computationally ex-pensive and not suitable for repeated runs necessary for a design and optimization process, also interpretation of the results can be difficult.

Koiter [1] introduced a perturbation approach where a perturbation expansion of the load parameter is made around a known state (e.g. a bifurcation buckling point) to capture the essential nonlinear behavior of a structure in the following form:

λ = λc+ aλcξ + bλcξ2+ . . .

In the above equation, ξ is the perturbation parameter, λ is the load parameter, λc

is the bifurcation buckling load, the coefficient a is known as the postbuckling slope, and the coefficient b is known as the postbuckling curvature. These postbuckling coefficients (a and b) can predict whether the postbuckling response is stable or unstable and can give a direct measure of imperfection sensitivity (see Fig. 1.1). When several buckling modes are associated with the critical bifurcation point,

(26)

the ﬁnal outcome of Koiter’s perturbation approach is a reduced set of nonlinear algebraic equations describing the relation between load parameter and perturbation parameters, where perturbation parameters are the unknowns.

Adaptation of Koiter’s perturbation approach within a finite element framework can remedy the afore mentioned problems in standard nonlinear analysis and is the main focus of this research. The resulting reduced order model leads to a significant decrease in the number of degrees of freedom involved in the finite element calcu-lation. Additionally, the approach can provide more insight into the fundamental characteristics of the structural behavior and can help in interpreting the results obtained by full model nonlinear analysis. Compared to standard full model non-linear analysis, such an approach is less accurate but still can be very effective in the preliminary design phase or in the optimization process where repeated analysis runs are necessary, but the required level of accuracy can be relatively lower.

(a) a6= 0 (unstable) (b) a = 0, b > 0 (sta-ble)

(c) a = 0, b < 0 (un-stable)

Figure 1.1: Typical postbuckling responses. The solid and dotted lines indicate, respectively, the response of the perfect and imperfect structures.

The present research belongs to a more general research line that aims at the development of fast tools for multi-fidelity nonlinear finite element analysis of struc-tures. The term ‘multi-fidelity’ refers to a hierarchical approach where analysis is carried out in a step by step manner increasing the level of analysis complexity at each step. The fast tools in this research line stand at an intermediate level of com-plexity between the linear analysis and the fully nonlinear analysis and bridge the gap between linear and full model nonlinear analysis as shown in Fig. 1.2. In special cases, semi-analytical or semi-empirical reduced complexity methods are available to bridge the gap. Semi-analytical methods are very fast and suitable for paramet-ric studies, but limited to problems of specific geometry and boundary conditions. Hence a more general finite element based reduced complexity method is desired and is attempted to achieve in the current research, which is a continuation of the previous research done by Tiso [2] in the same research group. It should be noted

(27)

1.2. LITERATURE REVIEW 3 that, as shown in Fig. 1.2b, the ﬁnite element based reduced complexity method is not proposed as a substitute for full model analysis. Rather they are meant to complement each other.

(a) Current approach in nonlinear analysis (b) Proposed hierarchical, multi-ﬁdelity approach for nonlinear ﬁnite element analysis

Figure 1.2: Comparison between the current approach and the proposed approach.

In recent years efforts have been made by several research groups to apply per-turbation methods in finite element context, which will be discussed along with some involved specific issues in the following section. Based on that discussion, in the subsequent section the research objectives of the present research will be formulated.

1.2 Literature Review

In the following a brief review of the relevant literature concerning the current re-search will be presented. Therefore, in light of the discussion made in the preceding section, perturbation methods with emphasis on ﬁnite element implementation will be the main focus of the discussion. The review is further divided into three parts based on the speciﬁc topics covered in the current research, namely single mode initial postbuckling analysis, multi mode initial postbuckling analysis and extension to dynamics, where dynamic buckling and nonlinear vibration problems are empha-sized. Although not treated in the present research, perturbation methods can also be applied for general transient problems [3–7].

(28)

1.2.1 Single mode initial postbuckling analysis

In this section, literature review basically on the finite element adaptation of Koi-ter’s perturbation approach in a single mode context will be discussed. The approach itself will be discussed in detail in the next chapter. Koiter-type perturbation ap-proach has been traditionally used in analytical and semi-analytical context [8–14]. In the context of finite element formulation an early research and review in this field was done by Gallagher [15] who outlined some difficulties involved in the finite element implementation of Koiter’s perturbation approach. They include, firstly, a lack in demonstration of the feasibility in large scale, complex practical problems and secondly, accurate computation of higher order modes and postbuckling coeffi-cients. Later research however addressed these problems, which is discussed in the following.

As already mentioned, accurate computation of b coeﬃcient is a key issue. Pig-nataro et al. [16] showed the diﬃculties involved in the accurate computation of

b coeﬃcients. It was found that geometrically exact nonlinear beam models are

necessary in order to obtain a correct b coefficient, and most of the technical beam theories where certain kinematical terms are neglected may lead to qualitatively wrong estimate of b coefficients. Salerno and Lanzo [17], and Pacoste and Eriks-son [18] proposed a geometrically exact beam based on Antman’s [19] kinematical model and obtained correct b coefficient. In a later work Pacoste and Eriksson [20] and recently Garcea [21] used beam elements based on corotational formulation to solve the problem. Lanzo [22] proposed a kinematic simplification to Green-Lagrange strain measure known as ‘Simplified Lagrangian’ that partly alleviates the problem. For instance, in case of Euler’s strut problem it gives zero b coefficient, which is of course not correct but the effect of the error remains globally small in the sense that, it results in an increment of 1% of load for a midspan displacement of 10% of the length. Concerning this issue, Lanzo [22] also argued that in case of thin walled plate and shell type structures due to their internal hyperstaticity, buckling is followed by a strong redistribution of the in-plane stress and geometrically exact kinematics are not necessary to obtain correct b coefficient.

Another important issue is the locking problem in terms of convergence of b co-efficient. The cause of the problem lies in the different degree of approximation of in-plane and lateral displacements which often holds for many beam and shell ele-ments causing inaccurate prediction of postbuckling stresses and consequently it also affects b coefficients. Different researchers adopted different techniques to circum-vent the problem. Two major contributions in this area were made by Olesen and Byskov [23, 24], and Poulsen and Damkilde [25]. Olesen and Byskov used Lagrange multipliers to introduce additional fields in the formulation. Poulsen and Damkilde solved the problem by including a local stress contribution leading to correct evalua-tion of postbuckling stresses. In order to tackle the problem, Casciaro [26] employed bubble function, and Lanzo et al. [27] used high continuity element, while Salerno and Lanzo [17], Pacoste and Eriksson [18, 20] and Tiso [2] used strain measures which are constant over the element.

(29)

1.2. LITERATURE REVIEW 5 important issue. This kind of locking affects the correct evaluation of buckling load (in case of linear buckling analysis carried out from undeformed configuration) and postbuckling coefficients which is dominant when the structure is characterized by high axial/flexural stiffness ratio in the presence of moderate prebuckling rotations. Garcea et al. [28] showed that perturbation algorithm based on independent ex-trapolation of displacement and stresses (i.e. mixed formulation) provides a good solution to this problem. In [28] a beam element was used with mixed formulation. Later Garcea [29] extended the formulation for high continuity shell element. Garcea et al. [28, 29] also showed that ‘frozen geometry’ approach in which prebuckling ro-tations are neglected may work well in many cases.

The perturbation approach used in this research is valid asymptotically in the neighborhood of the starting point of the perturbation expansion. All the works discussed above considered up to second order expansion of the displacement ﬁeld which is often enough to capture the initial postbuckling response. However in order to increase the range of validity of the perturbation expansion one can con-sider further higher order terms. Damil and Potier-Ferry [30], Azrar et al. [31, 32] and Lopez [33, 34] considered such higher order terms and showed the increase of the range of validity of the perturbation expansion for plate and shell problems. Following this approach Cochelin [35] proposed a path following technique where perturbation approach is applied in a stepwise manner i.e. the perturbation expan-sion is made at several points in the nonlinear solution path. Cochelin also proposed the possibility of using Noor’s [36] reduction method within the same framework. However in multi mode context considering such higher order terms will result into computation of too many postbuckling coeﬃcients which may not be suitable from the viewpoint of computational cost.

The literature discussed until now considered problems characterized by bifurca-tion type buckling. Initially Haftka et al. [37] and later Carnoy [38] and Casciaro et al. [39] considered snap through buckling problem where buckling occurs at a limit-point. The approach assumes the existence of a fictitious perfect structure having a bifurcation point and the behavior of the real structure is predicted by considering an ‘intrinsic’ imperfection in the fictitious perfect structure. The approach worked well with some beam structures but not with the same level of success in case of shell structures. Carnoy [38] proposed several measures for possible improvements. Sig-nificant contributors in the finite element implementation of perturbation approach, have already been mentioned (within a single mode context). Few more works that are worth mentioning are Eckstein et al. [40], Wu and Wang [41], Geier [42].

1.2.2 Multi mode initial postbuckling analysis

Multi mode study is necessary to assess the degrading eﬀect of load-carrying capac-ity due to modal interaction among clustered buckling modes. Within an analytical framework, modal interaction issue for discrete model was extensively treated by Thompson and Hunt [43], while Koiter [44, 45] considered it for continua. Following the same research line two more signiﬁcant analytical works were done by Srid-haran [46] and Kasagi [47]. Detailed derivation of the perturbation approach for

(30)

multi mode case can be found in the works of Byskov and Hutchinson [48], Peek and Kheyrkhahan [49], Salerno and Casciaro [50]. The derivation of Salerno and Casciaro is valid for simultaneous or nearly simultaneous buckling modes, while that of Byskov and Hutchinson, and Peek and Kheyrkhahan are valid for simultaneous, nearly simultaneous or well separated modes. Peek and Kheyrkhahan considered prebuckling nonlinearity and in contrast to others, made the expansion from any arbitrary point of the fundamental path (although the ﬁrst bifurcation point was preferred).

Using the expansion [50] several modal interaction studies were done by Lanzo and Garcea [22], Casciaro et al. [51], Salerno and Casciaro [50] and Garcea [29]. They used used high continuity ﬁnite elements [27, 52] and considered folded plate structures, for instance, T, C, L, Box cross-section beams (transversely loaded) and columns (axially loaded). Salerno and Casciaro [50] addressed issues like secondary bifurcations, post-critical attractive paths and worst imperfection shape [53, 54]. Garcea et al. [55] showed that the secondary bifurcation paths obtained using the perturbation approach are also in good agreement (in the neighborhood of the bifur-cation point) with standard full model path following approach, provided the same bifurcation branch is followed in the full model analysis. In the same line of research Bilotta et al. [56] considered shear deformability and composite plates, while Salerno and Lanzo [17] considered plane frames using an appropriate beam element. Modal interaction in folded plate type structures was also extensively studied by Menken et al. [57,58], Erp and Menken [59,60] and Schreppers and Menken [61]. They used the expansion given by Byskov and Hutchinson [48] and used a spline ﬁnite-strip type element. Based on the expansion given in [49] Kheyrkhahan and Peek [62] studied modal interaction in isotropic cylindrical and spherical shells. A review of works on modal interactions in thin-walled structures can be found in [63].

The multi mode postbuckling analyses discussed so far are basically reduced complexity analysis in the initial postbuckling regime. However this approach can also be used in conjunction with full model path following analysis, where the per-turbation approach carried out at any bifurcation point gives a predictor for the path following through any speciﬁc bifurcating branch. In other words perturbation approach is used for branch switching purpose. Some of the interesting works in this ﬁeld have been done by Kouhia and Mikkola [64], Huang and Atluri [65], and Magnusson [66].

1.2.3 Extension to dynamics

In the present research, the extension of the formulation to dynamics includes dy-namic buckling and nonlinear vibration problems. Therefore, the discussion of this section is also confined to dynamic buckling and nonlinear vibration problems. As far as dynamic buckling is concerned, according to Budiansky and Roth [67], the dynamic buckling load is defined as the level at which a large increase occurs in the deflection amplitude when the nonlinear equations of motion of the system are solved for different load levels. Recently Mallon [68] considered the dynamic buckling problem of a base-exited thin beam and cylindrical shell with top mass using a

(31)

semi-1.3. RESEARCH OBJECTIVES AND SPECIFIC ISSUES 7 analytical approach along with experimental validation. Some other analytical and semi-analytical works in this field include Petry and Fahlbusch [69], Aksogan and Sofiyev [70], and Jabareen and Sheinman [71]; while some finite element based works include Saigal et al. [72] and Yaffe and Abramovich [73]. In the context of reduced complexity analysis, Koiter’s perturbation approach for static postbuckling analy-sis can also be used to study dynamic buckling problems by taking inertial effects into account. Budiansky [74] devised the analytical foundation for this extension. Following the approach proposed by Budiansky, Schokker et al. [75] investigated dynamic instability of interior ring stiffened composite cylindrical shells under hy-drostatic pressure using a p-version of finite element in a single mode context. Using the same approach in a multi mode context, recently Chen and Virgin [76] considered various static and dynamic aspects of post-buckled thin plates both under mechan-ical and thermal loading. They used the same high continuity element [52] that was used previously by Lanzo and Garcea [22] for their static perturbation analysis.

Perturbation approach can also be used effectively for nonlinear vibration prob-lems. Under large excitations, a structure can exhibit a complicated dynamic be-havior. When finite displacements occur, it is appropriate to include geometrically nonlinear effects which cause the dependence of vibration frequency on the vibra-tion amplitude. The nonlinear vibravibra-tion behavior of structures, and of shells in particular, is therefore a topic that receives continuing attention [77, 78]. Many studies in the field of nonlinear vibration make use of a semi-analytical (Galerkin or Rayleigh-Ritz) formulation in order to describe the spatial dependence of the solution, however they are restricted to structures with a relatively simple geom-etry [79]. Rehfield [80] applied Koiter type perturbation approach for nonlinear (large amplitude) undamped free vibration problem in an analytical context and derived analytical expression for dynamic b coefficient (bD) for hinged straight beam

and simply supported rectangular plates. Wedel-Heinen [81] applied perturbation approach to study the influence of an initial geometrical imperfection on the un-damped, free, linear vibration (small amplitude) frequency of a structure at a given conservative load. This study was also done in an analytical context and simply supported beam and plate problems were treated. In the present research group, Jansen [82, 83] combined and extended the work of Rehfield and Wedel-Heinen and applied to vibration analysis of anisotropic cylindrical shell in a semi-analytical context, while Hermens [84] implemented Wedel-Heinen’s approach and Tiso [85] implemented Rehfield’s approach in a finite element framework. Azrar et al. [32] also applied perturbation approach for large amplitude free vibration problem of thin elastic plates in finite element context. An additional contribution of this work is the inclusion of up to 18th order terms in the perturbation expansion in order to increase the range of validity.

1.3 Research Objectives And Speciﬁc Issues

The aim of the present research is to show the feasibility of implementing the per-turbation approach in a general purpose ﬁnite element code in order to handle

(32)

prac-tical shell structures. In accordance with this perspective, the implementation of the present work will be done in the development environment of a general purpose finite element code DIANA [86] so that practical problems can be addressed. DIANA is a rare general purpose finite element code where a finite element implementation of Koiter’s initial postbuckling theory is available. The present implementation makes use of DIANA’s original implementation and does the necessary extensions. The following research objectives are set:

• To carry out primarily, single mode (only one buckling mode is considered,

usu-ally the first one) static initial postbuckling analysis of general shell structures. Accurate determination of postbuckling coefficients is emphasized. A key as-pect is the inclusion of prebuckling nonlinearity in the finite element based perturbation approach, which makes it possible to consider an important class of problem such as axially loaded cylindrical and conical shells, where buckling and postbuckling behavior are strongly influenced by the prebuckling nonlin-earity. Another key aspect is the incorporation of the capability to analyze composite shells, which are commonly used in practical structures.

• To extend the single mode analysis to multi mode (several buckling modes are

considered, usually the clustered modes) to account for the degrading eﬀect of modal interactions. The problem of selecting appropriate modes will be addressed.

• To make further extension to account for inertial eﬀect so that dynamic

buck-ling and vibration problems can be considered.

• To demonstrate the applicability and eﬃciency of the approach in reference

and bench mark problems. The possibility of application of the perturbation approach in the optimization process of variable stiﬀness composite panels will be explored.

Several specific issues involved in the finite element adaptation of perturbation methods have been discussed in the literature review section. In the following, stand point of the present research with respect to these issues is stated briefly. In the current research, two DIANA elements, namely a Mindlin type three-node, curved beam element and a eight-node, curved shell element are used. Regarding exact kinematics, in the present research ‘Simplified Lagrangian’ (see Section 1.2.1) approach is adopted for the Mindlin type three-node, curved beam element. In case of plate and shell type structures where due to their internal hyperstaticity, buckling is followed by a strong redistribution of stress, Lanzo [22] argued that, geometrically exact kinematics is not required to obtain correct b coefficients. Following this argument, the existing kinematics (see Eq. (A.2) in Appendix A) is used for the shell elements.

Unlike Kirchoff type elements none of the DIANA elements used in this research shows the locking effect in terms of convergence of b coefficient. In these elements the axial/in-plane and lateral/out-of-plane displacements are interpolated by quadratic

(33)

1.4. THESIS LAYOUT 9 shape functions. Possibly, because there is no difference in the degree of the interpo-lation polynomial, the locking problem related to b coefficient does not occur. Apart from these two elements, the shell element used by Tiso [2] is enhanced for compos-ite materials and implemented in DIANA, with the same constant strain measure over the element and it has been possible to obtain correct b coefficient without any locking problem.

Concerning extrapolation locking, in the present research, when prebuckling non-linearity is not too high ‘frozen geometry’ approach is adopted and when prebuck-ling nonlinearity is significant, prebuckprebuck-ling path is followed by a full model nonlinear analysis and perturbation analysis is carried out from the deformed configuration which is close to the bifurcation buckling load and the effect of prebuckling nonlin-earity is considered in the way proposed by Cohen [10], Fitch [11] and Arbocz and Hol [87].

For multi mode case, the perturbation expansion given by Byskov and Hutchin-son is used because it is valid for simultaneous, nearly simultaneous and well sep-arated buckling modes. Both for single mode and multi mode cases the eﬀect of prebuckling nonlinearity is included.

1.4 Thesis Layout

The thesis layout is briefly described in this section. Following this section, in Chapter 2, single mode initial postbuckling analysis is considered. In this chapter the perturbation approach and its finite element implementation is discussed in a single mode context with emphasis on the inclusion of prebuckling nonlinearity. The special case of linear prebuckling state is also discussed. Both Mindlin and Kirchoff type shell elements are used. Postbuckling coefficients are compared between these two type of elements, with semi-analytical results [14] and results available in the literature. In Chapter 3, the single mode formulation of Chapter 2 is extended to a multi mode context where expansion proposed by Byskov and Hutchinson [48] is used and the additional terms resulting from prebuckling nonlinearity are included. Sev-eral benchmark/reference problems made of both isotropic and composite material are considered. Issues including modal interaction, selection of appropriate modes, worst imperfection shapes are addressed and results are compared with full model nonlinear static analysis. The efficiency of the reduced order model is also demon-strated by comparing the computational time with full model analysis. Chapter 4 is oriented towards a specific applications. In this chapter basically the potential of the perturbation approach in the optimization process for the improvement of postbuckling performance of variable stiffness composite panels is demonstrated. In Chapter 5, the formulation presented in Chapters 2 and 3 for initial postbuckling analysis is extended with the inclusion of inertial effects and dynamic buckling and nonlinear vibration problems are addressed. Finally, in Chapter 6, conclusions of the present research is drawn and some recommendations for future work in this research line are proposed.

(34)

(35)

Chapter 2 Single mode initial postbuckling

analysis

2.1 Introduction

This chapter deals with a finite element formulation of Koiter’s perturbation ap-proach for single mode initial postbuckling analysis. In this apap-proach a perturbation expansion for the load and the displacement field is made around the bifurcation buckling point. The first order terms in the perturbation expansion of the displace-ment field are indeed the bifurcation buckling modes, while the second order terms are known as second order modes. Often a limited number of buckling modes and second order modes can be used to describe the initial postbuckling behavior. The number of equations in the resulting reduced set of nonlinear algebraic equations is the same as the number of buckling modes chosen in the perturbation expan-sion, and the perturbation approach can be regarded as the basis of a nonlinear reduced order model. Assuming that a single mode is associated with the lowest bifurcation buckling load, together with the computation of the buckling and second order mode, the postbuckling slope (a coefficient) and curvature (b coefficient) are computed. These postbuckling coefficients directly give a measure of the stability and imperfection sensitivity of the structure. For instance, in the case of conical and cylindrical shells one has a zero a coefficient and typically a negative b coef-ficient indicating unstable postbuckling behavior and imperfection sensitivity. The buckling mode, second order mode, and the corresponding postbuckling coefficients are characteristics of the perfect structure. Once they are computed the effect of geometric imperfections can be evaluated with very little additional computational cost for various imperfection shapes and magnitudes using the reduced order model, which in the single mode case consists of one nonlinear algebraic equation.

As discussed in the literature review section of the first chapter, often prebuckling state is assumed to be linear in the finite element formulation of Koiter’s pertur-bation approach. For an important class of problems consisting of cylindrical and conical shells under axial compression the assumption of a linear prebuckling state often leads to an overestimation of the bifurcation buckling load, and the buckling mode can be strongly different from the one obtained including a nonlinear

(36)

ling state. Therefore, the eﬀect of prebuckling nonlinearity will be included in the present study.

Many of the earlier works on Koiter’s theory are based on the principle of sta-tionary potential energy. However, in this thesis an alternative procedure will be used proposed by Budiansky and Hutchinson [88], that expresses the field equations directly in variational form using the principle of virtual work. In Budiansky and Hutchinson’s work the prebuckling state was assumed to be linear. Cohen [10] and Fitch [11], and later Arbocz and Hol [14, 87] derived the modifications necessary to include prebuckling nonlinearity. The present work makes use of the derivations done by Arbocz and Hol [14, 87] within a finite element context.

Two different types of shell elements are used. The first type consists of a class of Mindlin type (considers transverse shear stress) curved shell elements available in DIANA. Out of that class of elements, the 8 node, quadrilateral type element named CQ40S and CQ40L (for isotropic and laminated composite material, respectively, see Appendix A) are used in the numerical examples presented in this chapter. The another type of element that is used is a Kirchhoff [89, 90] type (transverse shear stress is ignored) triangular flat shell element developed by Allman [91–93]. Tiso [2, 94] enhanced this element to alleviate locking problem concerning conver-gence of b coefficient. In the present work the same element is further enhanced for composite material and implemented in DIANA element T18SH (see Appendix B). Past research showed that for accurate computation of the b coefficient, particularly with Bernoulli type beam and Kirchhoff type shell elements special element formu-lation is required [22, 23, 25, 27, 51, 94]. It is an additional finding that an accurate computation of the b coefficient is possible without any modification in the element formulation in case of Mindlin/Reissner [90, 95, 96] type beam and shell elements used in this work. In order to demonstrate this, a curved 3D beam element, CL18B available in DIANA, is used in one of the numerical examples.

In this chapter problems with both linear and nonlinear prebuckling behavior are treated and both isotropic and composite plate and shell structures are consid-ered. The postbuckling coefficients will be compared using the semi-analytical tools ANILISA [14] and BAAC [97], and also with the reported postbuckling coefficients by other researchers. The efficiency of the reduced order model that is obtained is also illustrated. The initial postbuckling response and limit-point buckling load are compared with full model finite element analysis carried out using DIANA with the same element.

2.2 The Perturbation Method

In this section the perturbation method for buckling and postbuckling analysis will be discussed with the inclusion of prebuckling nonlinearity. A detailed derivation of the equations is available in the report by Arbocz and Hol [87]. Here the basic procedure will be explained mentioning the essential equations. The functional no-tation introduced by Budiansky [74] will be used. First, the perturbation approach for the perfect structure will be explained. Next, it will be extended for the

(37)

imper-2.2. THE PERTURBATION METHOD 13 fect structure. In the following, symbols with bold font denote vector and tensor quantities while the scalar symbols are written in normal font.

2.2.1 Functional notation

If u and ² denote displacement and strain ﬁelds respectively, then in functional notation strain-displacement relation is written as

² = L1(u) +

1

2L2(u) (2.1)

Here L1 is a linear operator and L2 is a quadratic operator. Therefore, L1(u) is

a linear functional representing the linear part of strain and L2(u) is a quadratic

functional representing the nonlinear part of strain. Further, the bilinear operator

L11 is deﬁned such that

L2(u + v) = L2(u) + 2L11(u, v) + L2(v) (2.2)

From Eq. (2.2) it follows that

L11(u, v) = L11(v, u) (2.3)

L11(u, u) = L2(u) (2.4)

In this work Green-Lagrange strain has been used and in Eq. (A.2) of Appendix A, explicit deﬁnition of the operators L1, L2 and L11 has been shown with respect to

this speciﬁc strain-displacement relation.

2.2.2 Perfect structure

Let u, ², f, and σ be the generalized displacement, strain, load, and stress variables. Then the nonlinear strain-displacement relation is given by Eq. (2.1) and the linear elastic constitutive relation can be written as

σ = H(²) (2.5)

where H is a linear operator. The equilibrium equation in variational form is written as

σ· δ² − f · δu = 0 (2.6)

Here σ· δ² and f · δu denote, respectively, the internal virtual work of the stress σ through the strain variation δ², and the external virtual work of the load f through the displacement variation δu, both integrated over the entire structure. Further, it follows from Eqs. (2.1) and (2.2) that the ﬁrst order strain variation δ² produced by δu can be written as

δ² = L1(δu) + L11(u, δu) (2.7)

Additionally, in case of linear elasticity the following reciprocity relation holds

(38)

In this study proportional loading is considered, i.e. f =λf0. Now the variables

(u, ², σ) of the postbuckling equilibrium state can be expanded in the following perturbation expansion series about the prebuckling equilibrium state (u0, ²0, σ0)

at the same value of the variable load parameter λ

u = u0(λ) + u1ξ + u2ξ2+ u3ξ3+ . . .

² = ²0(λ) + ²1ξ + ²2ξ2+ ²3ξ3+ . . .

σ = σ0(λ) + σ1ξ + σ2ξ2+ σ3ξ3+ . . .

(2.9)

The variables (u0, ²0, σ0) are assumed to be nonlinear functions of λ = λ(ξ), while

the expansion functions (uk, ²k, σk) where k = 1, 2, . . ., are independent of λ and

ξ. The perturbation expansions (2.9) are assumed to be asymptotically valid in the

neighborhood of the bifurcation point deﬁned by λ = λc and ξ = 0.

Substituting Eqs. (2.9) into Eqs. (2.1), (2.5), and (2.6), taking the limit ξ → 0 and with some further manipulations one obtains the necessary equations for the bifurcation buckling load λc and the corresponding buckling mode u1

²1 = L1(u1) + L11(uc, u1) (2.10)

σ1 = H(²1) (2.11)

σ1· δ²c + σc· L11(u1, δu) = 0 (2.12)

where the subscript ()c denotes prebuckling quantities evaluated at λ = λc.

Next, it is assumed that the prebuckling variables can be expanded in the Taylor series u0 = uc+ (λ− λc) ˙uc + 1 2(λ− λc) 2_u_¨ c + . . . σ0 = σc+ (λ− λc) ˙σc+ 1 2(λ− λc) 2_σ_¨ c + . . . (2.13) where ˙( ) = _∂λ∂ ( ).

In addition it will be assumed that (λ− λc) admits the asymptotic perturbation

expansion

λ− λc = aλcξ + bλcξ2+ . . . (2.14)

In view of Eq. (2.14), if a plot of the load parameter λ versus the mode amplitude

ξ is made then a and b coeﬃcients respectively indicate the slope and curvature of

the postbuckling curve (see Fig. 1.1 of Chapter 1). In the present study symmetric bifurcation with zero postbuckling slope, a = 0, and typically negative postbuckling curvature, b < 0, indicating unstable postbuckling behavior is considered.

Inserting equations (2.13) and (2.14) together with equation (2.9) into equations (2.1), (2.5) and (2.6) and equating the coeﬃcients of ξ2 _{with the assumption of a = 0}

one can ﬁnally obtain the necessary equations for the determination of the second order mode u2

²2 = L1(u2) + L11(uc, u2) +

1

2L2(u1) (2.15)

(39)

2.2. THE PERTURBATION METHOD 15

σ2· δ²c+ σc· L11(u2, δu) + σ1· L11(u1, δu) = 0 (2.17)

In order to obtain the expression for the b coeﬃcient one can set δu = u1 in

Eqs. (2.12) and (2.17) and can make use of the reciprocity relation (2.8). This gives

b =−(1/λc∆)ˆ {2σ1· L11(u1, u2) + σ2· L2(u1)} (2.18)

where

ˆ

∆ = 2σ1· L11( ˙uc, u1) + ˙σc· L2(u1) (2.19)

In case a6= 0, Eq. (2.18) takes a more complicated form

b =−(1/λc∆)ˆ {2σ1· L11(u1, u2) + σ2· L2(u1)

+ a ˆ∆[ ˙σc· L11(u1, u2) + σ1· L11( ˙uc, u2) + σ2· L11( ˙uc, u1)]

+ (1/2)(a ˆ∆)2[2σ1· L11(¨uc, u1) + ¨σc· L2(u1)]}

(2.20)

According to Cohen [10] Eq. (2.20) can be simpliﬁed by making u2 orthogonal

to u1 by applying the following orthogonality condition

˙

σc· L11(u1, u2) + σ1· L11( ˙uc, u2) + σ2· L11( ˙uc, u1) = 0 (2.21)

Therefore, Eq. (2.20) simpliﬁes to

b =−(1/λc∆)ˆ {2σ1· L11(u1, u2) + σ2· L2(u1)

+ (1/2)(a ˆ∆)2[2σ1· L11(¨uc, u1) + ¨σc· L2(u1)]}

(2.22)

2.2.3 Imperfect structure

Now it will be shown how the behavior of the imperfect structure can be derived from the properties of the perfect structure. If the initial geometric imperfection is denoted by ¯ξ ˆu where ¯ξ is the imperfection amplitude and ˆu is any arbitrary geometric

imperfection pattern, then the strain-displacement Eq. (2.1) can be modiﬁed to

² = L1(u) +

1

2L2(u) + ¯ξL11(ˆu, u) (2.23) Further, the asymptotic expansion as deﬁned by Eq. (2.14) is also modiﬁed to

ξ(λ− λc) = aλcξ2+ bλcξ3− αλcξ¯− β(λ − λc) ¯ξ + . . . (2.24)

where the coeﬃcients α and β are known as the ﬁrst and second imperfection form factors respectively. This equation and Eq. (2.14) can be seen as the reduced order model for the imperfect and perfect structure respectively. Using the same approach as explained in Section 2.2.2 the expressions for α and β are obtained as

α = (1/λc∆)[σˆ 1· L11(ˆu, uc) + σc· L11(ˆu, u1)] (2.25) β =(1/ ˆ∆){σ1 · L11(ˆu, ˙uc) + ˙σc· L11(ˆu, u1) + H[L11( ˙uc, u1)]· L11(ˆu, uc) − αλc £ σ1· L11(¨uc, u1) + (1/2) ¨σc· L11(u1, u1) + H[L11( ˙uc, u1)]· L11( ˙uc, u1) ¤ } (2.26)

(40)

where ˆ∆ is deﬁned by Eq. (2.19).

Now, with the background established so far in this section one can do the postbuckling analysis by computing up to the second order terms in the perturbation expansion of the displacement ﬁeld u in Eq. (2.9). It can be achieved through the following step by step procedure:

1. Computation of the prebuckling state uc at λ = λc.

2. Computation of the buckling load λc and the corresponding buckling mode u1.

3. Computation of the second order mode u2.

4. Computation of the a and b coeﬃcients.

5. Computation of the ﬁrst and second imperfection form factors α and β. 6. Computation of the mode amplitude ξ corresponding to the applied load

pa-rameter λ using Eq. (2.24).

7. Recovering the displacement ﬁeld by substituting the already computed terms in Eq. (2.9).

2.3 Finite Element Implementation

2.3.1 Nonlinear buckling analysis

According to Eq. (2.7) the strain variation at the critical point (²c) can be written

as

δ²c = L1(δu) + L11(uc, δu) (2.27)

Insertion of Eqs. (2.10) and (2.11) into Eq. (2.12) together with Eq. (2.27) gives

H[L1(u1) + L11(uc, u1)]· [L1(δu)

+ L11(uc, δu)] + σc· L11(u1, δu) = 0

(2.28)

After some algebraic manipulation of Eq. (2.28) and replacing the L1 and L11

oper-ators and the continuous displacement ﬁelds u1, uc, δu respectively, with the ﬁnite

element matrices BL, BN L and nodal displacements q1, qc, δq one can have the

discretized form of Eq. (2.28)

δqT[BT_LHBLq1+ BTN L(qc)HBLq1+ BTLHBN L(qc)q1

+ BT_{N L}(q_c)HBN L(qc)q1+ B

T

N L(q1)σc] = 0

(2.29)

Because δq is an arbitrary displacement vector, remaining quantities inside the bracket are zero

BT_LHBLq1+ B T N L(qc)HBLq1+ B T LHBN L(qc)q1 + BT_{N L}(q_c)HBN L(qc)q1+ B T N L(q1)σc = 0 (2.30)

A perturbation approach for geometrically nonlinear structural analysis using a general purpose finite element code