Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 41, 3, pp. 561-574, Warsaw 2003 GEOMETRICALLY NONLINEAR STATIC ANALYSIS OF SANDWICH PLATES AND SHELLS Jakub Marcinowski Institute of Structural Engineering, University of Zielona Góra e-mail: J.Marcinowski@ib.uz.zgora.pl Sandwich shells composedof three layers: two thin andvery strong faces, and the soft and comparatively weaker core which is much thicker than their neighbouring layers covering it from up and down, are considered in the paper. The proposal of a finite element adequate to such a kind of structural members is presented. In the paper the finite element based principally on the Ahmad original element and on the author’s adap- tation of this very element to the nonlinear range is presented. It was assumed that in the faces and in the core thematerials exhibit orthotro- pic properties, and only elastic deformations are taken into account.The static analysis is performed within fully geometrically nonlinear range. The ultimate purpose of the work is determination of the critical value of the load intensity factor in a quasi static stability analysis. All proce- dures related to tracing of nonlinear equilibriumpaths are adopted from the codes prepared before for homogeneous shells. Two verification pro- blems are inserted. They confirm correctness of the adopted approach. Keywords: sandwich shells, othotropy, nonlinear stability analysis, finite element method 1. Introduction Sandwich shells andplates exhibitmanyvaluable properties and this is the reason that structures made in that technology are so often used in engine- ering practice. Themost important advantages of sandwich structures can be listed as follows: good thermal and acoustic isolation, good vibration damping, good strength toweight ratio, good local and global buckling resistance.These 562 J.Marcinowski physical andmechanical properties of sandwich structural elements have cau- sed that they are applied not only as finishing elements but also as structural members (cf. Hop, 1980; Romanów, 1995). The mechanical behaviour of this kind of structures is analysed in this work, and particular emphasis is put on static analysis with reference to the equilibriumstability phenomenonwithin the range of largedisplacements.The main purpose of the present work can be explained with the help of Fig.1. The critical value of the load intensity factor for one parameter loading is searched. This particular value can be found in fully nonlinear, quasi static analysis in which nonlinear equilibrium paths are determined in the whole load displacement space. The critical points in the form of limit points or bifurcation points can appear on these paths, and the primary critical point defines the searched value of load parameter. The procedure of calculation of the nonlinear equilibrium paths must be comprehensive enough to determine all mentioned objects no matter how complicated the primary or secondary paths could be. Fig. 1. On parameter loading and equilibrium path Theproblem is important and this is the reason that somany authors have been trying to build better and better numerical tools to describe themecha- nical behaviour of these kinds of structural members. In the literature there exist many proposals of numerical modelling of mechanics of sandwich shells by FEM. In this context, it is worthy to mention the work of Rammerstorfer et al. (1992). The proposed models are more and more comprehensive and take into account phenomena observed in experiments. Ferreira et al. (2000) present the approach in which the plasticity is taken into account and degene- rated conception of finite element is used. Riks and Rankin (1995) proposed Geometrically nonlinear static analysis... 563 a manner of description of face delamination. Works of Ding and Hou (1995) andMoita et al. (1999) refer to initial buckling of sandwich structures. In the present work the conception of a degenerated finite element is used to model the mechanical behaviour of sandwich, three layered shells within the range of assumptions formulated in the next section. 2. Main assumptions A three layered sandwich shell or plate is considered. The curvature of the shell can be arbitrary. The structure is composed of a thick and soft core co- vered from both sides by very thin and strong faces (comp. Fig.1a).Materials of both are linearly elastic and exhibit orthotropy properties. The directions of orthotropy can be different in the core and in faces. Themain assumptions refer to the mode of core deformation. It has been assumed that the straight normal to the middle surface of the core remains straight but can rotate independently with respect to the deformed middle surface. It means that the severe Kirchhoff’s assumption is dropped. The faces can deform only within their planes and these membrane defor- mations are equal to their counterparts in the core fulfilling in this way the strain continuity conditions. Local deformations of the faces, independent of the mode of deformation of the core are excluded from considerations, which is the obvious drawback of the presented approach. Themodes of deformations taken into account are shown in Fig.2a,b. The mode shown inFig.2c is an example of themode excluded fromconsiderations. Fig. 2. Modes of deformation Generally, arbitrary displacements are taken into account but as far as rotations are concerned they should be moderate in the presented version of the program.Strains are small, so linear stress-strain relations can be adopted. The strain energy corresponding to the stress normal to the middle surface is ignored. It means that the plane stress state is assumed within every layer parallel to the middle surface. The constant, plane, membrane stress state occurs within both faces. 564 J.Marcinowski The formulatedassumptions correspondto theassumptionsof theMindlin- Reissner theory of shells extended to the case of sandwich shells. 3. Finite element In the present approach the original conception of degeneration due to Ahmad et al. (1970) is adopted. The idea is based on such geometry and displacement approximation that the reduction of 3D solid to the 2D curved surface fulfils all assumptions of theMindlin-Reissner shell theory. Fig. 3. Geometry approximation The geometry approximation is defined as follows (see Fig.3)    x y z    = ∑ i Ni(ξ,η) 1+ ζ 2    xi yi zi    t + ∑ i Ni(ξ,η) 1− ζ 2    xi yi zi    b (3.1) Geometrically nonlinear static analysis... 565 where: {x,y,z}⊤ are global coordinates of any point within the shell element volume, {xi,yi,zi} ⊤ t(b) – global coordinates of the top (t) or bottom (b) node, ξ,η,ζ – curvilinear coordinates of any point within the shell element volume, Ni(ξ,η) – the shape function corresponding to the node i, and summation holds over all nodes of the element. Relation (3.1) can be rewritten in the following alternative form    x y z    = ∑ i Ni(ξ,η)    xi yi zi    m + ∑ i Ni(ξ,η)    (V3i)x (V3i)y (V3i)z    (3.2) where {xi,yi,zi} ⊤ m are global coordinates of the node lying on the middle surface, {(V3i)x,(V3i)y,(V3i)z} – components of the nodal vector. It is worthy to mention that the parent element for the defined one is the cube of side 2. It means that {ξ,η,ζ} ∈ 〈−1,1〉. It is apparent from relations (3.1) and (3.2) that the defined approximation is linear with respect to the third coordinate ζ. The above definitions refer to the whole volume: the core together with two faces. Fig. 4. Displacement approximation The displacement approximation is more general, and is defined as follows    u v w    = ∑ i Ni(ξ,η)    ui vi wi    m + ∑ i Ni(ξ,η)ζ ti 2    (υ̂1i)x −(υ̂2i)x (υ̂1i)y −(υ̂2i)y (υ̂1i)z −(υ̂2i)z    { αi βi } (3.3) where {ui,vi,wi} ⊤ are components of nodal displacements in the global co- ordinate system, {u,v,w}⊤ are components of displacements at a given point {ξ,η,ζ}, [υ̂1i,−υ̂2i] – matrix composed of unit vectors defining axes of the nodal coordinate system (cf. Fig.4), {αi,βi} ⊤ – two independent rotations defined appropriately in the nodal coordinate system. 566 J.Marcinowski Fig. 5. Local coordinate system In this definition formula, the linear approximation of displacements with respect to the third coordinate is visible. It is the guarantee that the plane section will remain plane after deformation. Two independent rotation para- metersmakepossible free rotation of the straight normal to themiddle surface, according to the adopted assumptions. In this manner, the displacement field within the whole sandwich element is described by means of five nodal parameters: three translations and two independent rotations. It is exactly the sameas itwasdone inoriginalAhmad’s element. To define strains, a local coordinate system (LCS) is introduced in such a way that the axes x′ and y′ occur in the plane parallel to themiddle surface, and the axis z′ is perpendicular to it. Strain-displacements relations ensue from the full Green-de Saint Venant strain tensor and take the following form    εx′ εy′ γx′y′ γx′z′ γy′z′    =    ∂u′ ∂x′ + 1 2 [(∂u′ ∂x′ )2 + (∂v′ ∂x′ )2 + (∂w′ ∂x′ )2] ∂v′ ∂y′ + 1 2 [(∂u′ ∂y′ )2 + (∂v′ ∂y′ )2 + (∂w′ ∂y′ )2] ∂u′ ∂y′ + ∂v′ ∂x′ + ∂u′ ∂x′ ∂u′ ∂y′ + ∂v′ ∂x′ ∂v′ ∂y′ + ∂w′ ∂x′ ∂w′ ∂y′ ∂u′ ∂z′ + ∂w′ ∂x′ + ∂u′ ∂x′ ∂u′ ∂z′ + ∂v′ ∂x′ ∂v′ ∂z′ + ∂w′ ∂x′ ∂w′ ∂z′ ∂v′ ∂z′ + ∂w′ ∂y′ + ∂u′ ∂y′ ∂u′ ∂z′ + ∂v′ ∂y′ ∂v′ ∂z′ + ∂w′ ∂y′ ∂w′ ∂z′    (3.4) The component εz′ was omitted due to the assumption that the influence of σz′ on the strain energy is ignored. Only first three components will appear within the faces of the sandwich element. Geometrically nonlinear static analysis... 567 The constitutive relations are formulated in LCS as well, and due to the assumed orthotropy properties, take the form    σx σy τxy    =    Ex 1−νxyνyx Exνxy 1−νxyνyx 0 Eyνyx 1−νxyνyx Ey 1−νxyνyx 0 0 0 Gxy       εx εy γxy    (3.5) within the faces, and    σx σy τxy τxz τyz    =    Ex 1−νxyνyx Exνxy 1−νxyνyx 0 0 0 Eyνyx 1−νxyνyx Ey 1−νxyνyx 0 0 0 0 0 Gxy 0 0 0 0 0 kGxz 0 0 0 0 0 kGyz       εx εy γxy γxz γyz    (3.6) inside the core. In these relations the following denote: Ex,Ey – Young’s moduli in the x and y directions, νxy – ratio of the lateral strain |εx| to the strain εy when loaded in the y direction, Gxy,Gxz,Gyz – shearmoduli in the planes xy, xz and yz respectively; the condition Exνxy =Eyνyx must hold. The correction factor k = 5/6 was introduced here due to the fact that the uniform shear stress distribution follows from relation (3.3) while it is parabolic actually. The governing equations of the problem are obtained from the principle of virtual work. If p denotes the vector of external load this principlewritten for the whole sandwich shell takes the following form ∫ VC (σC)⊤δεC dv+ ∫ VTF (σTF)⊤δεTF dv+ ∫ VBF (σBF)⊤δεBF dv− ∫ A p ⊤δu dA=0 (3.7) where VC, VTF ,VBF are volumes of the core, of the top face and the bottom face respectively, A is the area where the external load is applied, δε denotes strain variations due to virtual displacements, δu are virtual displacements. Passing to finite elements, this principle takes the following form ∑ (e) ( ∫ VCe σiδεi dv+ ∫ ATF σiδεi dA+ ∫ ABF σiδεi dA−λFiδdi ) =0 (3.8) 568 J.Marcinowski where: tt, tb are thicknesses of the top and bottom faces respectively, λ – load intensity factor, Fi – nodal forces due to the load on the reference level, V Ce – volume of the core, ATF ,ABF – area of the top and bottom face respectively,∑ (e) – denotes aggregation over all elements, di are nodal displacements. Twotermscorresponding to the strain energygenerated in the faces appear in this relation. Because the strains are definedby the same nodal parameters, the above relation can be rewritten, after some further derivations, as follows ∑ (e) [ (KC +KTF +KBF)de−λFe ] =0 (3.9) where: KC,KTF ,KBF are stiffness matrices of the core, top face and bottom face, respectively, Fe – vector of the nodal forces due to the load on the reference level. All stiffness matrices are dependent on the nodal displacements. The resulting set of nonlinear algebraic equations assumes the following form Ψ(d,λ)=KNd−λF =0 (3.10) where KN is the global stiffness matrix dependent on d, d – global vector of the nodal displacements, F – global vector of the nodal forces. It is worthy tomention that KC is calulated using the reduced integration scheme 2×2×2 Gauss points, while KTF , KBF are calculated using 2×2 integration scheme (2D integrals). To obtain a load-displacement curve, set (3.10) have to be solved for the whole range of the load intensity factor λ. The adopted procedures are exactly the same like these described in thework byMarcinowski (1999) and need not to be discussed here. These procedures are versatile enough to trace even the most complicated equilibrium paths with bifurcation points, limit points, secondary paths, etc. 4. Verification problems To verify the correctness of the proposed approach and the computer pro- gram itself, a few problems were solved. The first one was taken from the library of verification problems of the COSMOS/M system (cf. [2]). It refers Geometrically nonlinear static analysis... 569 to deflections of a square sandwich plate clamped along all edges and loaded by a uniformpressure. The geometry andmaterial data are given in Fig.6, on which the load intensity factor versus central deflection curve was shown. The presented solution was comparedwith that obtained byCOSMOS/Mand the solution taken from the paper of Schmit and Monforton (1970). Quite good correspondence can be observed. The core in this problem exhibits only shear rigidity, while the faces are rigid in the plane and in the lateral direction. As a matter of fact Gxz and Gyz of the faces could not be incorporated in the present model because only the membrane effect has been taken into account within the faces. Fig. 6. Nonlinear equilibrium path for the square clamped plate As a second example the well known benchmark of Sabir and Lock (1972) will be considered. It is a cylindrical shell simply supported along their rec- tilinear edges and loaded laterally by a concentrated force applied to its cen- ter. Details of boundary conditions, geometry and material data are shown in Fig.7. In this figure FE mesh is also shown. This division is chosen after checking that the shell does not exhibit bifurcation phenomenon. 570 J.Marcinowski Fig. 7. Cylindrical shell loaded by a concentrated force Fig. 8. Equilibrium paths for the homogeneous shell Geometrically nonlinear static analysis... 571 The homogeneous case of such a shell was solved by many authors. A comparison of the present solution with those obtained by other authors is shown in Fig.8. Very good correspondence can be observed. Using the same geometry, a new problem is created following the idea of Ferreira et al. (2000). The thickness of the shell is now divided into three zones. Two thin faces (tf = 0.635mm) are separated from the homogeneous shell and, in this way, a sandwich shell is created. To analyse the influence of the core rigidity on the overall stiffness of the sandwich shell the Young modulus of the core (and the bulk modulus which is expressed by E, and G = E/[2(1+ ν)]) was reduced 10, 100 and 1000 times, while the material parameters of the faces are kept constant. In this manner, three cases of a more andmoreweak core are considered. The results of analysis are presented in Fig.9 (load versus vertical displacements of the point C) and in Fig.10 (load versus vertical displacements of the point A). They are shown together with the solutions obtained bymeans of the systemCOSMOS/M for the same data. Fig. 9. Equilibrium path for the sandwich shell Discrepancies can be observed for very large displacements. These displa- cements are accompanied byfinite rotationswhichwere not taken into account 572 J.Marcinowski Fig. 10. Equilibrium path for the sandwich shell in the presented version of program. It is the reason of these discrepancies. The coincidence is pretty good within the initial portions of the curves, and particularly in the vicinity of the critical points. 5. Final remarks The conception of degeneration usually adopted to homogeneous shells has turned out to be effective also in the case of sandwich shells. The perfor- med test confirmed that the proposed approach is correct within the adopted assumptions. The limit ratio tf/tc (the face thickness over the core thick- ness) was established, and it seems that the adopted assumptions are valid for tf/tc < 1/15. There is possibility of extending the proposed approach to a case of thicker layers and amultilayer case, but it would require integration along the lateral direction within every layer. There is nopossibility of taking into account local deformations of the faces (cf. Fig.1c), and this is the obvious drawback of the presented approach. All numerical procedures used before for nonlinear static analysis of homogeneous plates and shells (cf. Marcinowski, 1999) turned out to be effective also in the case of sandwich shells. Geometrically nonlinear static analysis... 573 References 1. Ahmad S., Irons B.M., Zienkiewicz O.C., 1970,Analysis of thick and thin shell structures by curvedfinite elements, Int. J. Num.Meth. Engg., 2, 419-451 2. COSMOS/M2.5, Electronic documentation, Structural Research andAnalysis Corporation, Los Angeles 1999 3. Crisfield M.A.,1981, A fast incremental/iterative solution procedures that handles ”snap-through”,Computer and Structures, 13, 55-62 4. Ding Y, Hou J., 1995, General buckling analysis of sandwich constructions, Computer and Structures, 55, 485-493 5. FerreiraA.J.M., Barbosa J.T.,MarquesA.T., De Sa J.C., 2000,Non- linear analysis of sandwich shells: the effect of core plasticity, Computer and Structures, 76, 337-346 6. Hop T., 1980,Konstrukcje warstwowe, Warszawa, Arkady (in Polish) 7. Marcinowski J., 1999, Nieliniowa stateczność powłok sprężystych, Wrocław, OficynaWydawn. PolitechnikiWrocławskiej (in Polish) 8. Moita J.S., Soares C.M.M., Soares C.A.M., 1999, Buckling and dyna- mic behaviour of laminated composite structures using a discrete higher-order displacementmodel,Computer and Structures, 73, 407-423 9. Rammerstorfer F.G., Dorniger K., Starlinger A., 1992, Composite and sandwich shells, In: Nonlinear analysis of shells by finite elements, Edit. F.G. Rammerstorfer, Springer-Verlag,Wien-NewYork, 131-194 10. Riks E., Rankin Ch.C., 2002, Sandwich modelling with an application to the residual strength analysis of a damaged compression panel, International Journal of Non-linear Mechanics, 37, 897-908 11. Romanów F., 1995, Wytrzymałość konstrukcji warstwowych, Zielona Góra, Wydawn.Wyższej Szkoły Inżynierskiej, (in Polish) 12. SchmitL.A.,MonfortonG.R., 1970,Finite deflectiondiscrete element ana- lysis of sandwich plates and cylindrical shells with laminated faces, AIAA Jo- urnal, 8, 8/9, 1454-1461 13. Sabir A.N., Lock A.C., 1972, The application of finite elements to the large deflectiongeometricallynonlinearbehaviourof cylindrical shells, In:Variational Methods in Engineering, Sauthampton Univ. Press, 7166-7175 14. Stein E., BergA.,WagnerW., 1982,Different levels of nonlinear shell the- ory in finite element stability analysis, In: Buckling of Shells, Edit. E. Ramm, Proceedings of a State-of-the-Art Colloquium, Springer-Verlag, Berlin Heidel- berg NewYork 574 J.Marcinowski Geometrycznie nieliniowa analiza statyczna płyt i powłok warstwowych Streszczenie W pracy rozważa się płyty i powłoki warstwowe złożone z trzech warstw: dwóch twardych okładek i stosunkowo miękkiego rdzenia, który jest znacznie grubszy od przykrywającychgo okładek.Zaprezentowanoelement skończony, opartynakoncepcji degeneracji ośrodka trójwymiarowego, dobrze opisującywłasnościmechaniczne takiej powłoki. W opracowanym elemencie wykorzystano oryginalną koncepcję Ahmada, Ironsa i Zienkiewicza rozszerzoną przez autora na zakres geometrycznie nieliniowy. Zakładając błonowy stan odkształceń i naprężeń w okładkach, uzupełniono macierz sztywności rdzenia o macierze sztywności okładek wykorzystując przy tym te same stopnie swobody, które wprowadził Ahmad. Założono, żemateriały rdzenia i okładek są liniowo sprężyste i wykazują własności ortotropii. Przedstawiona analiza obejmuje zagadnienia quasistatyczne w zakresie dużych przemieszczeń i umiarkowanych obro- tów ze szczególnym uwzględnieniem utraty stateczności równowagi. Do wyznaczania geometrycznie nieliniowych ścieżek równowagiwykorzystanoprocedury autorskie sto- sowanewcześniej z powodzeniemwprzypadkupowłokpełnych, jednorodnych.Wpra- cy przedstawiono kilka przykładówpotwierdzających poprawność proponowanej kon- cepcji. Przeprowadzone testy wykazały, że sformułowane założenia i oparta na nich koncepcja opisu są prawdziwe dla powłokwarstwowych,w których stosunek grubości okładki do grubości rdzenia jest mniejszy od 1/15. Manuscript received December 2, 2002; accepted for print March 11, 2003