70_

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 2, 419-451 (1970) ANALYSIS OF THICK AND THIN SHELL STRUCTURES BY CURVED FINITE ELEMENTS SOHRABUDDIN AHMAD* Civil Engineering Department, University of Engineering and Technology, Dacca, East Pakistan BRUCE M. IRONS? AND 0. C. ZIENKIEWICZ$ Civil Engineering Department, University of Wales, Swansea SUMMARY A general formulation for the curved, arbitrary shape of thick shell finite elements is presented in this paper along with a simplified form for axisymmetric situations. A number of examples ranging from thin to thick shell applications are given, which include a cooling tower, water tanks, an idealized arch dam and an actual arch dam with deformable foundation. A new process using curved, thick shell finite elements is developed overcoming the previous approximations to the geometry of the structure and the neglect of shear deformation. A general formulation for a curved, arbitrary shape of shell is developed as well as a simplified form suitable for axisymmetric situations. Several illustrated examples ranging from thin to thick shell applications are given to assess the accuracy of solution attainable. These examples include a cooling tower; tanks, and an idealized dam for which many alternative solutions were used. The usefulness of the development in the context of arch dams, where a ‘thick shell’ situation exists, leads in practice to a fuller discussion of problems of foundation deformation, etc., so that practical application becomes possible and economical. INTRODUCTION The analysis of shells with an arbitrarily defined shape presents intractable analytical problems. If, in addition, the shell is a thick one in which shear deformation is significant the applicability of classical approaches become questionable. In civil engineering structures ranging from water tanks and cooling towers to pressure vessels and arch dams, these difficulties must be overcome if satisfactory (and hopefully optimized) designs are ever to be achieved. Two paths are at present possible-structural models and numerical analysis by computer. Both present their own difficulties but clearly for design purposes the latter is more convenient. In formulating numerical solutions to elastic shell problems many alternative approximations may be used. That of the finite element approach will be followed here. In idealization of the shell by finite elements a geometrical simplification of replacing the curved shell by an assembly of flat elements is most frequently ~ s e d . l - ~ With this simplification a large number of elements must invariably be used, and any advantage that can be gained by more sophisticated elements (which in smaller numbers can yield improved accuracy) is lost. Thus the need for elements which can take up curved shapes becomes obvious in this context. * Formerly Commonwealth Research Scholar, Civil Engineering Department, University of Wales, Swansea. t Lecturer. $ Professor of Civil Engineering and Head of the Department. Received 1 October 1969 419 @ 1970 by John Wiley & Sons, Ltd. 420 SOHRABUDDIN AHMAD, BRUCE M. IRONS AND 0. C. ZIENKIEWICZ Some attempts to develop such curved shell elements have been made in recent years, but to date their application is limited to shallow shell situations and to those in which shear deformation is The development of large curvilinear elements for three-dimensional analysis by the relatively simple process of ‘isoparametric’ formulation (coupled with numerical integration) appears to open a possible avenue: Elements of the type shown in Figure 1 have for some years been used Figure 1. Three-dimensional hexahedral elements of parabolic and cubic types with success for three-dimensional analysis purposes.*-ll As the dimensions of such elements are completely arbitrary (and indeed are obtained simply by specifying nodal co-ordinates) one could visualize their ‘attenuation’ to represent geometrically a prescribed shell segment. Indeed such a representation would forcefully bring it home to the engineer that the classification of problems into shells, plates, etc. is an artificial one-introduced merely for the sake of con- venience. Obviously in nature only three-dimensional problems exist and the sub-groups are introduced to reduce analytical or computational labour. With a straightforward use of the three-dimensional concept, however, certain difficulties will be encountered. In the first place the retention of three degrees of freedom at each node leads to large stiffness coefficients for relative displacements along an edge corresponding to the shell thickness, This presents numerical problems and inevitably leads to ill-conditioned equations when shell thick- nesses become small compared with the other dimensions in the element. The second factor is that of economy. The use of several nodes across the shell thickness ignores the well-known fact that even for thick shells the ‘normals’ to the middle surface remain practically straight after deformation. Thus an unnecessarily high number of degrees of freedom has to be carried, involving penalties of computer time. In this paper a specialized formulation is presented overcoming both these difficulties.12 The constraint of straight ‘normals’ is introduced to overcome the second problem and the strain THICK AND THIN SHELL STRUCTURE 42 1 energy corresponding to stresses perpendicular to the middle surface is ignored. With these modifications an efficient tool for analysing curved thick shells becomes available. Its accuracy and wide range of applicability is demonstrated on several examples. The reader will note that the two constraints introduced correspond only to a part of the usual assumptions of shell theory. Thus, the statement that after deformation the normals remain normal to the deformed middle surface has been deliberately omitted. This omission permits the shell to experience shear deformations-an important feature in thick shell situations. Further, as integration is carried out numerically it is not necessary to introduce the various simplifying assumptions always present in conventional shell theory and resulting in the familiar situation of a wide variety of differential equations, apparently different, but describing the same problem. Indeed, in the very nature of the process, the need for formulating such equations disappears. The paper is divided into three parts. The first is devoted to the description of the basic theory. While at this stage it is unnecessary to repeat the fundamentals of the ‘displacement’ approach to finite element theory, adequately given in a text,’ it is felt essential to put on record some of the special mathematical manipulations needed to put the previously described assumptions into practice. Tt is hoped that this section will permit the suitably equipped reader to repeat and implement the computer programs used. While the first part contains the ‘meat’ of the presentation the reader may proceed on cursory reading directly to the second and third parts presenting applicatidns of the analysis. The second is concerned mainly with verijication and outlines several examples of varying complexity for which solutions are available. The third deals in some detail with the special problems of arch dams-or indeed other similar shell structures in which foundation interaction is present. THEORY Geometric definition of the element Consider the two typical thick shell elements of Figure 2(a). The external faces of the element are curved, while the sections across the thickness are generated by straight lines. Pairs of points, itop and ibottom, each with given Cartesian co-ordinates, prescribe the shape of the element. Let c, 7 be two curvilinear co-ordinates in the middle plane of the shell and 5 a linear co- ordinate in the thickness direction. If further we assume that 5, 7, 5 vary between - 1 and I on the respective faces of the element we can write a relationship between the Cartesian co-ordinates of any point of the shell and the curvilinear co-ordinates in the form Here Ni(.$,q) is a function taking a value of unity at the node i and zero at all other nodes. If the basic functions Ni are derived as ‘shape functions’ of a ‘parent’ element, square (or triangular) in plan (Figure 2b) and are so ‘designed’ that compatibility is achieved at interfaces, then the curved space elements will fit into each other. It is well known that for the first kind of element the functions are parabolic and for the second cubic in 5 and q and thus the curved shape of the shell element can take up a parabolic or cubic form respectively. By placing a larger number of nodes on the surfaces of the element more elaborate shapes can be achieved if so desired. Suitable shape functions for the elements of Figure 2(b) are listed in Appendix I. The relation between the Cartesian and curvilinear co-ordinates is now established and it will be found desirable to operate with the curvilinear co-ordinates as the basis. 422 SOHRABUDDIN AHMAD, BRUCE M. lRONS AND 0. C. ZIENKIEWICZ It should be noted that the co-ordinate direction 5 is only approximately normal to the middle surface. Figure 2. General curved shell elements; (a) Parabolic and cubic thick shell elements; (b) Parabolic and cubic ‘parent’ elements ; (c) Geometry, local co-ordinates and nodal displacements It is convenient to rewrite relationship (1) in a form specified by the ‘vector’ connecting the upper and lower points (i.e. a vector of length equal to the shell thickness t ) and the mid-surface co-ordinates. Thus we have (Figure 2c) with THICK AND THIN SHELL STRUCTURE 423 Displacement jield The displacement field has now to be specified for the element. As the strains in the direction normal to the mid-surface will be assumed to be negligible, the displacement throughout the element will be taken to be uniquely defined by the three Cartesian components of the mid-surface node displacement i and two rotations of the nodal vector V,$ about orthogonal directions normal to it. If two such orthogonal directions are given by vectors fzi and 8, (of unit magnitude) with corresponding (scalar) rotations ai and pi we can write, dropping the suffix ‘mid’ of equation (2): + where u, u and w are displacements in the directions of the global x , y and z axes. As an infinity of vector directions normal to a given direction can be generated, a particular scheme has been devised to ensure a unique definition. This is given in Appendix 11. Once again if Ni are compatible functions then displacement compatibility is maintained between adjacent elements. It can also be shown that the displacement definition reproduces any state of rigid body motion-a condition necessary for convergence.*l? lo Physically, it has been assumed in the definition of equation (3) that no strains occur in the direction 5 . While this is not exactly normal to the middle surface it represents to a good approximation one of the usual shell assumptions. At each node i of Figure 2(c) we have now the five basic degrees of freedom. Dejnition of strains and stresses To derive the basic properties of a finite element the essential strains and stresses have to be defined.’ The components in directions of orthogonal axes related to the surface 5 = constant are essential if account is to be taken of the basic shell assumptions. Thus if at a point in this surface we erect a normal z’ with two other orthogonal axes x’ and y’ tangent to it (Figure 2c) the strain components of interest are ( E ’ } = (4) with the strain in direction z’ neglected so as to be consistent with the shell assumption. It must be noted that in general none of these directions coincide with those of the curvilinear co- ordinates t, 7, t, although x‘, y’ are in the t - ~ plane (5 = constant). * As the definition of the co-ordinates is more general than that of the displacements but includes it as a special case the elements are called superparametric, to distinguish from the isoparametric elements of References 8-1 I . 424 SOHRABUDDIN AHMAD, BRUCE M. IRONS AND 0. C. ZIENKIEWICZ E [D’] = - (1 -9) The stresses corresponding to these strains are defined by a matrix (Q’} and are related by the elasticity matrix [D’]. Thus - l v 0 0 0 1 0 0 0 - 0 0 1 - v V - 0 I-v 2k 1-v 2k - sym. where {E;} may represent any ‘initial’ strains due, for instance, to thermal expansion. The 5 x 5 matrix [D’] can now include any anisotropic properties and indeed may be prescribed as a function of 5 if sandwich construction is used. For the present moment we shall define it in which E and v are Young’s modulus and Poisson’s ratio respectively. The factor k included in the last two shear terms is taken as 1-2 and its purpose is to improve the shear displacement approximation. From the displacement definition it will be seen that the shear distribution is approximately constant through the thickness, whereas in reality the shear distribution is approximately parabolic. The value k = 1-2 is the ratio of relevant strain energies. Element properties and necessary transformations over the volume of the element, which are quite generally of the form The stiffness matrix-and indeed all other ‘element’ property matrices-involves integrals where the matrix [S] is a function of the co-ordinates. In the stiffness matrix [SI = [BIT PI PI (4 = [Bl{V for instance, with the definition (7) so that [B] relates the strains to the nodal parameters. The theory of the subject is given in the text’ and many other references, and need not be dwelt upon here. If in the present context we can express [S] as an explicit function of the curvilinear co-ordinates and transform similarly the infinitesimal volume, dx dy dz, then a straightforward (numerical) integration will allow the properties to be evaluated. Thus some transformations are necessary. Equation (3) relates the global displacements u, u, w to the curvilinear co-ordinates. THICK AND THIN SHELL STRUCTURE 425 The derivatives of these displacements with respect to the global x , y , z co-ordinates are given by a matrix relation [J] = au av aw I - ay & 5 ax ay aZ a7 5 Zj ax ay aZ - a { Z - - [ - au - av - aw aZ ax aZ The Jacobian matrix is defined as = [JI-' and is calculated from the co-ordinate definition of equation (2). Now, for every set of curvilinear co-ordinates the global displacement derivatives caii be obtained numerically. A further transformation to local displacement directions x', y' , z' will allow the strains, and hence the [B] matrix, to be evaluated. First the directions of the local axes have to be established. A vector normal to the surface 5 = constant can be found as a vector product of any two vectors tangent to the surface. Thus Following the process which defines uniquely two perpendicular vectors, as given in Appendix 11, and reducing to unit magnitudes, we construct a matrix of unit vectors in x', y', z' directions (which is in fact the direction cosine matrix) [dl = [%,%%I (1 3) The global derivatives of displacements u, u and w are now transformed to the local derivatives of the local orthogonal displacements by a standard operation au aU aw 426 SOHRARUDDIN AHMAD, BRUCE M. IKONS AND 0. C. ZIENKIEWICZ From this the components of the [B‘] matrix can now be found explicitly The infinitesimal volume is given in terms of the curvilinear co-ordinates as dx dy dz = determinant [J] d ( d7 d 5 (16) and this standard expression completes the basic formulation. The computer programs use Gaussian quadrature for the integration. The two-point rule suffices in the 5 direction, while a minimum of three or four points in both 6 and 71 directions is needed for parabolic or cubic elements respectively. Some remarks on solution The element properties are now defined, and the assembly and solution are standard processes. A particularly effective scheme of equation-solving known as the ‘front solution’ is used because it is inherently better than a banded solution for elements containing mid-side nodes.13 It remains to discuss the presentation of the stresses, and this problem is of some consequence. The strains being defined in local directions, {a’} is readily available. These are indeed directly of interest but as the directions of local axes are not easily visualized it is convenient to transfer the components to the global system using the following expression If the stresses are calculated at a nodal point where several elements meet then they are averaged. In a general shell structure, such as a doubly curved arch dam, the stresses in a global system do not, however, give a clear picture of shell surface stresses. The matrix-vector handling scheme14 (developed for convenience in dealing with the various transformations-see Appendix 111) therefore includes a ‘Jacob? eigenvalue routine which diagonalizes the stress tensor giving the principal stresses. The direction cosines of such stresses are obtained as vectors. Regarding the shell surface stresses more rationally, one may note that the shear components rX,:.,, and T~,: . , , are in fact zero there, and can indeed be made zero at the stage before converting to global components. The values directly obtained for these shear components are the average values across the section. The maximum transverse shear value occurs on the neutral axis and is equal to 1.5 times the average value. ‘VERIFICATION’ AND GENERAL EXAMPLES General remarks Most numerical processes are by their very nature approximate, because of curtailment of significant digits (round-off) or various physical idealizations introduced into the process of solution. An assessment of accuracy attainable with the present computational scheme, therefore, THICK AND THIN SHELL STRUCTURE 427 needs.to be made. In this section a series of examples ranging from thin to thick shell situations will be presented to demonstrate the accuracy attainable and hence the limitations of the solution. Several special classes of problems will be discussed. Axi-symmetric shell with axi-symmetric load I f the geometry'of the shell is that of a body of revolution and in addition the loading is symmetric about the axis, the general program can be simplified considerably. Only two co- ordinates are now necessary to define the geometry of a particular point and the degrees of freedom at a typical node reduce to three instead of five. Details of such modifications are fully described in Reference 12 and need not be repeated here. Several examples in this class of problem will be given. Example 1 Spherical dome under uniform pressure, Figure 3. An 'exact' solution based on thin shell theory15 is known for this simple case, In Figure 3 a comparison of moments and hoop forces is given for a subdivision using 24 cubic elements. The size of elements is graded, the smallest being used in the vicinity of the encastrk end. 01 I Figure 3. Spherical dome under uniform pressure analysed with 24 elements. (First element subtends an angle of 0.1" from the fixed end, others in arithmetic progression.) A44 = Meridional bending moment in.-lb/in.; T = Hoop force Ib/in.; v = Q 428 SOHRABUDDIN AHMAD, BRUCE M. lRONS AND 0. C. ZlENKIEWlCZ When the shell is thin, negligible shear effects are expected, but the differences obtained due to applying the pressure at inner and outer surfaces are worth noting. (The conventional thin shell theory applies to all loads at the middle surface.) The accuracy obtained is excellent. Example 2 Thin cylinder under radial edge load, Fig