J. Stegmann, J. C. Rauhe, L. Rosgaard and E. Lund Institute of Mechanical Engineering
Aalborg University, Denmark e–mail: js@ime.auc.dk
ABSTRACT
Summary This paper deals with a finite element formulation of a 4-node element for geometrically non-linear (GNL) analysis of layered shell structures. The element employs the MITC approach proposed by Dvorkin & Bathe [1] in order to eliminate the problem of shear locking. It uses full numerical integration both in-plane and through the thickness. The formulation passes all patch tests and shows good results in all tested configurations.
Introduction
The quest for robust finite elements for shell analysis has since the introduction of shell elements been a major goal. The goal seems to be almost reached with the introduction of the Assumed Natural Strain methods, which include the MITC (Mixed Interpolation of Tensorial Components) approach applied in the present work. A natural next step is to generalize these robust elements to allow the analysis of structures ranging from conventional metal structures, over composite laminates to sandwich structures. Furthermore, the use of these enhanced elements in other fields such as design optimization seems an obvious development. Such work has already begun and the present work is to a great extend related to the work of e.g. [2, 3, 4] but further introduces robust shell finite elements into the field of topology optimization. These robust shell element are essential for future design and optimization of advanced structures such as wind-mill wings.
Element formulation
Considering some of the typical applications of laminated shells, such as ship hulls, wind-mill wings, aeroplanes and other large structures it seems reasonable to limit the element modelling capabilities to gross response of structures. Consequently, an equivalent single layer model for describing the laminate behavior seems adequate. Furthermore, considering the material properties and layer stacking of such large structures, it is noted that these often change rapidly over the structural surface. Consequently, the mesh of finite elements in such a model must be relatively fine to model the changing material properties, and in such situations the benefit of a coarse mesh of higher order elements is lost. Thus, it seems a sound startegy to choose a linear element which is also less expensive in computational cost. However, the purely displacement-based 4-node shell element suffers heavily from shear locking and consequently, the MITC method is incorporated in the formulation thus eliminating these problems.
The current element formulation is based on the degenerated solid approach and the kinematic assumptions enforced classifies the formulation as a Mindlin-type element with the following kinematic interpolation:
in which ais the a’th node,ui are the displacements,Na is the interpolation, tis the thickness coordinate andhais the element thickness at nodea. The first term of (1) represents the in-plane contribution and the second term the out-of-plane contribution. The latter is described by the dis-placement of the reference surface normal0
V a
3i
but this can be expressed through the rotation about the two in-plane node directors0V1iaand0V2iawhich are defined as shown in Fig. 1.
1
Figure 1: The 4-node shell element (a) and definition of unit node director0Va3in nodea(b). The rotations and about0V2iaand0V1iarespectively express the displacement of the surface normal.
From the kinematic interpolation the element strain-displacement matrix,B, may be formulated directly, but to eliminate shear locking (emanating from spurious transverse shear energies) rows 5 and 6 are corrected according to the MITC method. The MITC method is founded in the obser-vation that the transverse shearing strains are computed correctly at the mid-sides of the 4-node element. The idea is to correct the B-matrix so that transverse shearing strains are computed in these four mid-side pointsA D(called tying points) in stead of in the Gauss-points. These cor-rect values are then interpolated across the element through a new set of interpolation functions:
~
are the new interpolation functions in r and scorresponding to theij’th strain and thek’th tying point. From the expression in (2) the assumed strains (AS) are expressed from the directly interpolated (DI) strains at the mid-side points. The new interpolations must naturally fulfill the relation:
so that thek’th interpolation function assumes the value1in thek’th tying point and the value0in all other tying points. The new interpolations are of the same order as the standard isoparametric interpolation functions used in the element. This leads to a corrected strain-displacement matrix and in turn to an element stiffness matrix, K, by expressing the constitutive matrix,C, with the following coefficients in the material coordinate system:
C
Note that the constitutive matrix resulting from the coefficients in (4) will have zero-values in the third row and column. This is done to enforce the Mindlin-assumption of negligible transverse effects. As shear correction factor we use5=6for all but sandwich structures.
The global system of equations is then assembled using the equilibrium of work, expressed for any element of volumeV and areaA:
Z where rdenotes a load vector and the superscripts b, c and s denote body-, concentrated and surface loads respectively. The stress is expressed as the Second Piola-Kirchhoff stress tensor to ensure work conjugacy with the Green-Lagrange strain tensor for finite strain. The global gov-erning equations are obtained as the sum of (5) over all elements in the discretized domain. The element integration is carried out by a full numerical Gauss integration. The in-plane integration is carried out as usual but the thickness-direction integration is carried out as piecewise integra-tion over all layers in the element. Consequently, the thickness-direcintegra-tion coordinate,t, is corrected using the transformation in (6) so thattl
2[ 1;1]in each layer: wherehlis the thickness of thel’th layer and the sum of the first term of the parenthesis constitutes the total thickness of the preceding layers. Inserting the modified coordinate in the Gauss integra-tion represents an integraintegra-tion by substituintegra-tion and consequently, the derivative of the coordinate transformation in (6) must be found and multiplied to the integrand.
Numerical examples
The element passes all patch tests performed in constant curvature, shear, twist and membrane as also noted by [1]. As benchmark examples we use the laminated cylindrical shell shown in Fig. 2. The geometry appears in various articles but the current example was adopted from [5]
who developed an enhanced full layerwise 3-node triangular element for geometrically non-linear analysis.
Figure 2: Cylindrical laminated shell used as numerical example.
The cylindrical shell is hinged at the straight sides and free at the two curved ends. The entire geometry is modelled with side length L = 508 mm, radius of curvature R = 2540 mm and angle2'=0:2rad. The lay-up is a +45/-45 (measured from the globalx-axis) laminate of equal thickness layers and a total thickness of h = 12:4mm. The material properties are stated in the material coordinate system (1-axis aligned with the globalx-axis) as:E1 = 3:2993106 MPa,
E
=0:25. The results from the current element formulation are presented in Fig. 3 and compared
0 0,2 0,4 0,6 0,8 1 1,2
0 2 4 6 8 10 12
Center point displacement mm
LoadkN To & Liu [5]
Current Linear solution
Figure 3: Results of linear and non-linear analysis of the cylindrical laminated shell in Fig. 2.
to the results obtained with a 88 mesh of triangular elements by [5].
As can be seen the results correlate very well although the current formulation exhibits slightly stiffer behavior. The current formulation also performs very well in linear analysis and numerous other non-linear examples.
Topology optimization of laminated shell structures has recently been implemented using the present element formulation and preliminary results show good correlation with those obtained by e.g. [6]. The development of optimization capabilities will continue in the near future, both for topology and shape design optimization.
REFERENCES
[1] Eduardo N. Dvorkin & Klaus-Jürgen Bathe. A continuum mechanics based four-node shell element for general non-linear analysis. J. Eng. Comp., 1, 77–88, (1984).
[2] G. Alfano, F. Auricchio, L. Rosati & E. Sacco. MITC finite elements for laminated composite plates. Int. J. Num. Met. Eng., 50, 707–738, (2001).
[3] S. Klinkel, F. Gruttmann & W. Wagner. A continuum based three-dimensional shell element for laminated structures. Computers & Structures, 71, 43–62, (1999).
[4] B. Brank & E. Carrera. Multilayered shell finite element with interlaminar continuous shear stresses: a refinement of the Mindlin-Reissner formulation. Int. J. Num. Met. Eng., 48, 843–
874, (2000).
[5] C. W. S. To & M. L. Liu. Geometrically nonlinear analysis of layerwise anisotropic shell struc-tures by hybrid strain based lower order elements. Finite Elements in Analysis and Design, 37, 1–34, (2001).
[6] S. J. Lee, J. E. Bae & E. Hinton. Shell topology optimization using the layered artificial material model. Int. J. Num. Met. Eng., 47, 848–867, (2000).