Meshing, Element Types, and Element Shape Functions

CHAPTER

2

Meshing, Element Types, and Element Shape Functions

King H. Yang Wayne State University, Detroit, Michigan, United States

2.1 STRUCTURE IDEALIZATION AND DISCRETIZATION Before analyzing the response of a structure subjected to loading, an FE model must be created. As demonstrated in Table 1.2, the computer software only has access to the nodal coordinates, element connectivity, boundary conditions, and loading conditions from the input data set. Based on these data, the software forms the element stiffness matrices [k], assembles them into the structure stiffness matrix [K], applies the boundary and loading conditions, and then uses the Gauss elimination or an equivalent method to find the nodal displacements. In other words, it is the user’s responsibility to discretize the structure, numerically idealize each member, prescribe boundary conditions, and apply loading conditions. In Example 1.2, the whole bridge structure is discretized into five truss members. Each truss member is idealized into a 2-node, 1D bar element. For each element, the forming nodes need to be arranged in a specific order, based on the way the element stiffness matrix is formulated. In FE terminology, the word mesh is defined as the collection of nodes (which contain information related to geometric locations) and elements (which prescribe the order of connectivity between nodes). After the mesh is created, boundary and loading conditions are applied at corresponding nodes. Finally, the Gauss elimination method is used to calculate nodal displacements. An example of the numerical discretization and idealization is shown in Fig. 2.1, where the cross-sectional view A-A of a long 3D dam is idealized as a 2D plane strain problem before the cross section is discretized into a finite number of triangular elements. Clearly, developing a 3D FE model of the entire dam is an option, if there are no concerns regarding costs associated with model development or computational resources. Even then, it would take a long time to simulate the response of a large model, which may not be desirable. Since the strength of a dam depends on the weakest section of the entire dam, it is reasonable and computationally more efficient to analyze only a representative 2D section of the entire 3D dam. Assuming that the z-axis is set to be along the axial direction of the dam, we can see from a transverse slice of the dam that force is applied to the xey plane. Because Basic Finite Element Method as Applied to Injury Biomechanics. http://dx.doi.org/10.1016/B978-0-12-809831-8.00002-7 Copyright © 2018 Elsevier Inc. All rights reserved.

51

52

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

FIGURE 2.1 Left: A sketch of a long 3D dam used for water storage. Inset: The cross-sectional view A-A, representing the weakest cross section of the dam, can be used to analyze the structure integrity of this dam under load. Right: This cross section is discretized into a finite number of triangular plane strain elements. Further discussion regarding the advantages and disadvantages of using triangular element type is presented in Section 2.3.

a dam is much thicker along the z-direction than along the cross-sectional directions, it is understandable that strain components involving z-axis are much smaller than those involving x- or y-axes. In other words, it is safe to assume that εzz ¼ gyz ¼ gzx ¼ 0 for a plane strain element. As a side note, a plane stress element (which represents a plate with a much smaller dimension in the z-direction as compared to those in the x- and y-directions) subjected to biaxial loading along the x- or y-axes would have nearly zero stress involving the z-axis. In other words, szz ¼ syz ¼ szx ¼ 0 are assumed for a plane stress element. The equations needed to represent plane stress and plane strain elements are provided in Section 1.2.3, Eqs. (1.22) and (1.23), respectively. Unlike the example problem shown above, all real-world engineering problems are 3D. Before modern, high-performance computers became readily available, numerous theories were developed and reported to reduce the total number of DOFs so problems could be solved with the relatively low-speed and low corememory computers of that time period. Examples aimed to reduce the total DOFs include the use of the axisymmetric solution principle, which is defined as solving a 2D symmetric problem when the structure can be formed by rotating a 1D line about a single point or solving a 3D problem by rotating a 2D edge about an axis of rotation. For example, a 2D circular disk can be simplified by rotating a straight line 360 degrees about the origin, and a 3D cylindrical structure with varying diameters over the vertical length can be simplified by rotating a 2D plate 360 degrees about the vertical axis (Fig. 2.2).

2.1 Structure Idealization and Discretization

FIGURE 2.2 An axisymmetric approach can be used to model a circular disk by rotating a line element 360 degrees about the center or to study a 3D vase by rotating a slice made of 2D elements 360 degrees about its axis.

True axisymmetric problems require all the geometry, material, and loading conditions to be symmetric about the same point or axis. As an example, the axisymmetric approach can be used to solve an aircraft turbine blade loading under normal operations. However, when attempting to solve a problem related to a couple of pieces of debris hitting one of the turbine blades, axisymmetric 2D method will not work because the load is asymmetric. Therefore, despite the savings in computational resources gained by reducing real-world 3D problems into axisymmetric 2D problems, finely meshed 3D models are preferred by modern FE model developers. Another commonly used method to reduce computational time is to capitalize on symmetric conditions. For example, if a symmetric plate is loaded under symmetric plane stress loading conditions, only a quarter of a plate is needed to solve the problem. Again, in order to use this type of simplification, there needs to be a perfectly symmetric condition. Considering that in most real-world problems loading is asymmetric, and modern computers are equipped with large amounts of core memory and large clusters of high-speed processors, saving computational resources becomes a less critical issue, unless an extremely large model is presented. Another advantage of using a full model instead of simplification based on symmetric conditions is the ability to detect artifacts. If a problem was truly symmetric, then the results would be symmetric. If asymmetric results were found, it would be obvious that there were numerical artifacts associated with poor meshing. Another consideration related to idealizing a structure includes omissions of some details (such as small holes or fillets) that are insignificant for the particular

53

54

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

FIGURE 2.3 A finer mesh size is used near the region of interest (ROI), located at the lower left corner, to calculate the stress increase due to concentrated loading (not shown) near the ROI. For the areas away from the ROI, the mesh size is much larger, because the stress variation would be much lower. This approach is advantageous when the computational resources are limited. The disadvantage of this approach is that a new model needs to be developed each time the ROI is changed.

problem. By skipping these details, an FE model can be created more quickly. Additionally, in the early days FE models were developed for specific types of engineering analyses. For example, at one time it was common practice to use small-sized elements near the areas of interest while coarse meshes were implemented for all other elements (e.g., Fig. 2.3). Again, problems that called for special approaches were mainly due to lack of computational resources. A model created with one specific analysis in mind could only be used for finding a solution for that particular purpose. A new model would need to be created if the finely meshed area were no longer of concern for a new analysis type. With the aid of meshing software, modern, general-purpose FE models tend to be created using a uniformly distributed mesh with a very fine mesh size, so one model is sufficient to study different boundary and loading conditions. Fig. 2.4 shows a general purpose, detailed FE car mesh developed at the National Crash Analysis Center for vehicle crashworthiness evaluations. The model can be used to simulate frontal, side, and rear-end impacts, which is quite an advantage over creating three separate FE models for the three different impact directions. The next two sections describe the fundamental nature of nodes and elements in more detail.

2.2 Node

FIGURE 2.4 A cut-away view of a Ford Taurus whole car-FE model developed at the National Crash Analysis Center (NCAC) of the George Washington University.

2.2 NODE The required information for a node, also known as a grid point, is the spatial location in a global coordinate system. In a global coordinate system, the location of a point is defined in a rectangular or spherical space. Commonly used coordinate systems include the Cartesian coordinate system, polar coordinate system, and spherical coordinate system. At each node, any constrained DOFs can be prescribed, and the responding DOFs due to applied forces and moments can be calculated. In theory, all nodes possess six DOFs, that is, three translations along the x-, y-, and z-axes and three rotations about the x-, y-, and z-axes. Typically, when these six DOFs are constrained, they are denoted as 1, 2, 3, 4, 5, and 6, respectively. For example, a 1, 2, 3, 5 constraint means there are no translations along any of the three axes (denoted as 1, 2, 3) and no rotation is allowed about the y-axis (denoted as 5). For some element types, certain DOFs are prohibited from happening, and these prohibitions need to be carefully considered. For example, 1D, 2-node truss, spring, and cable elements all allow three translational DOFs, but no rotational DOFs are permitted in 3D analysis. On the other hand, a 2-node beam element, which is geometrically idealized the same way as a 1D truss or spring element, admits all six DOFs at each node. Similarly, a 4-node membrane element contains no rotational DOFs (i.e., only three translational DOFs are allowed), but a 4-node plate element allows all three translational DOFs and two rotational DOFs about the two axes in the plane of the element (i.e., rotation about the axis perpendicular to the element surface is not calculated). Similar to a 4-node plate element, a 4-node shell element in certain FEA software packages does not allow in-plane rotation, which is also known as the drilling DOF. The reason for neglecting the drilling DOF is that there

55

56

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

can be incompatibilities with adjacent elements. However, the in-plane rotational DOF is permitted in some FEA software packages. Lastly, an 8-node hexahedral brick element allows only three translational DOFs. All forces that load the FE model need to be applied through the corresponding nodal translational DOFs, and moments need to be applied at the corresponding nodal rotational DOFs. The 8-node, hexahedral brick elements do not permit any rotational DOFs, and therefore application of moments to this element type has no effect. Similarly, any translational displacement constraint (e.g., zero nodal displacements in the x-, y- and/or z-axes) and any rotational constraint can only be applied at the node of the corresponding nodal DOF. If one particular element type allows no rotational DOFs, then applying any constraints to those DOFs will have no effect on the FEA results. For example, prohibiting zero translational and rotational movements of an 8-node hexahedral brick element will yield the same results as prohibiting zero translational DOFs alone. In some cases, the loading condition is described as a fixed amount of motion (such as compressing a tube for a fixed distance) or a constant velocity (such as a car-FE model hitting a rigid wall at a speed of 50 km/h), these prescribed motions need to be applied to the node(s), as well. For ease of use, some software packages allow users to apply distributed forces to a surface area. Internally, these packages still need to redistribute the surface loads into nodal loads. For FE solvers without such an option, users need to distribute the surface loads according to the corresponding element shape functions, to be described in Chapter 6.

2.3 ELEMENT Elements are the basic building blocks used to idealize structures. Because a building block can be made of any material types (e.g., metal, wood, bone), information related to its material properties must be contained within the element. Geometrically speaking, elements can be in the form of lines (1D trusses or beams), areas (2D membrane and plate), or solids (tetrahedral or brick). They are formed by connecting nodes that are arranged in a specific order, based on the mathematical relationships (shape functions) among these nodes. A wrong sequence provided to the input data deck will result in erroneous results. Unless stated otherwise, all nodes forming a 1D element are arranged from left to right, and all nodes in a 2D element are arranged in a counterclockwise order. For a 3D brick element representing an 8-node, 2-layer solid element, the nodes that form the lower layer of the element (P1 to P4) are arranged first, in a counterclockwise order, and the upper layer (P5 to P8) is subsequently arranged in the same counterclockwise manner. Note that the starting points for the bottom (P1) and top (P5) layers initiate at the same corner but are located at different heights along the z-axis. The set of element shape functions in the FE method is a collection of multipurpose interpolation functions needed for determining physical values of non-nodal points based on known nodal coordinates within an element. These functions can also be used to redistribute uniformly or nonuniformly distributed surface loads to

2.3 Element

nodal loads, and to find stress/strain contours. The order of a shape function depends on how many nodes there are in an element. For example, a linear (first order) interpolation is sufficient for a 2-node, 1D bar element, while a quadratic (second order) interpolation is needed for a 3-node, 1D bar element. Limited by the low computational power available for FEA in the early days, there was a tendency to use higher order interpolation element types in order to lower the number of elements, which in turn, reduced the total number of DOFs for obtaining solutions on these limited-power computers. These complicated or special higher order elements are more difficult to formulate than lower order elements, despite the fact that they are capable of satisfying a higher order of continuity. A simple continuity of the displacement field is called a C0 continuity. In most structural mechanics and injury biomechanics problems, strain is the most sought after response variable. Because strain can be calculated from the first derivatives of the displacement fields, a C0 continuity is sufficient for most structural mechanics and biomechanics problems. There are other problems which require continuity of the first derivatives of the displacement fields across the element edges. This type of continuity is known as C1 continuity. For example, some plate and shell elements based on KirchhoffeLove theories require C1 continuity in order to be valid. Fig. 2.5 graphically explains the essence of C0 and C1 continuity.

FIGURE 2.5 The column on the left depicts an element with C0 continuity, that is, the displacement or
 deflection field is continuous, but the velocity or slope field du dx is not. For an element with C1 continuity, both the displacement/deflection and velocity/slope fields are continuous, but the acceleration/curvature field is not.

57

58

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

2.3.1 SIMPLEST ELEMENT TYPES With ever-increasing computational power, modern FE models are built with extremely large numbers of elements to ensure close resemblance of the structures of interest. As the size of the elements becomes smaller and smaller, the difference in solutions between lower and higher order elements diminishes to a point of insignificance, in the majority of cases. As such, present-day FE modelers prefer the use of the simplest types of element with very fine mesh resolutions. Because higher order elements are not commonly used in contemporary FEA, theories behind the element types used for possessing a higher order of continuity will not be discussed in this book; only the simplest types of element will be discussed. At an early design stage, a coarse mesh model may be useful for capturing trends in the magnitudes of response variables that occur due to parameter changes. Simulations can be run quickly when using such a coarse mesh, and hence many parametric studies can be completed in a short period of time. A finer mesh model, which guarantees numerical convergence, is eventually needed to obtain a more accurate solution than that obtained from a coarse mesh model. In the FE method, the word “convergence” is related to how much discrepancy exists between the solutions calculated by the FE model compared to those obtained analytically or experimentally when the mesh is continuously refined. Without convergence, we cannot be confident that the FEA represents the real-world scenario. The source of discrepancy may come from lack of iteration convergence and mesh convergence. Regarding iteration convergence, some FE solutions are obtained through iterative procedures (e.g., the Jacobi method, to be described in Section 7.2). The iterative procedures will not stop until the error between two successive iterations is smaller than a preset value. Using such a method, the numerical and analytical/ experimental solutions become closer with a higher number of numerical iterations. If the numerical analysis is allowed to run long enough, the two solutions will be nearly identical. Hence, a lack of iteration convergence can be easily resolved by reducing the allowable error between two successive iterations. In terms of mesh convergence, all FE models involve some degree of simplification when creating a mesh to represent the structure of interest. Obviously, a finer mesh will result in a closer representation of the real-world structure to be modeled. In this case, it is understandable that the more elements are used to develop a model, the more accurate the results will be. After several mesh refinements, the mesh is considered “converged” if the differences in solutions between the current refinement and its predecessor are smaller than a preset value. In some occasions, continuing increase of the mesh density actually moves the solution away from convergence. This phenomenon usually indicates that the model is not properly defined. One possible reason may come from oversimplification of the model. As mentioned previously, omitting some fine details (such as not modeling small holes) is a common practice to reduce the total DOFs. In this case, the high stress concentration that occurs near the small hole cannot be exhibited unless the hole is

2.3 Element

–.15

0

.15

.3

.45

.6

.75

.9

1.05

FIGURE 2.6 The stress distribution for a linear elastic plate subjected to a concentrated load is twice as high for the model meshed with rectangular elements as compared to models meshed with triangular elements. Note that only eight divisions are set up to highlight the differences in stress distributions when different element types are used. As such, the peak element stresses discussed in the previous paragraph cannot be displayed in this figure. Those values are directly output from the FE solver. Even then readers can see that the “high stress” regions are different when different element types are selected to model the plate.

explicitly represented in the mesh. Other possibilities include the selection of different (and sometimes wrong) element types. For example, using a beam element type is more efficient than other element types when modeling a beam (See Section 2.5.3, Fig. 2.21, and Table 2.4). In addition to the requirement for convergence, FE models created using different element types or arranged in different manners could result in different solutions. Fig. 2.6 shows a rectangular plate loaded at the center of the left edge by a point (concentrated) load and fixed at the right edge of the plate. Even though the stress contours for all three models look highly similar, the magnitudes of the peak stresses are quite different. For the model developed with quadrilateral elements, the peak stress (a dimensionless magnitude of 3.9) is found to be more than twice that of the model built with right-angled triangle elements (a peak magnitude of 1.8) and equilateral triangle elements (a peak magnitude of 1.7). Here, no units (such as MPa) are provided for the stress, because the purpose of this exercise is to find relative values when the same material properties are used for all three models, without referencing any particular material.

2.3.2 1D ELEMENT TYPE The 1D element type includes, but is not limited to, a 2-node linear or a 3-node quadratic element. jThe 2-node elements, considered to be the simplest 1D element type (or pseudo 2D or 3D line elements when they are not in line with global axes), include the structure types of truss, spring, cable, beam, and frame. The only load a 2-node truss member can support is axial (i.e., no bending load is allowed). As such, no resistance can be generated from a truss member when a transverse force (vertical load) is applied. However, when a number of truss members are connected together by pin joints, the structure as a whole can resist applied forces

59

60

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

FIGURE 2.7 A taper truss structure (left) can be idealized as several truss elements (right), each with a constant cross-sectional area.

that are not aligned with the axial direction of the truss member. Because truss members are pinned together, these connecting members are free to rotate about the axes of the pins, and hence there is no resistance to moments. Additionally, this structure type needs to have a constant cross-sectional area, a constant elastic modulus, and can be oriented in any direction. For idealizing a truss member with varying cross-sectional areas, as shown in Fig. 2.7 (left), a series of elements, each with a different constant cross-sectional area, need to be used (Fig. 2.7 right). The number of elements used for idealizing a truss is based on the desired degree of accuracy. A spring can only support an axial load (tension or compression) and has only one DOF for each node along the axis of the spring. Sometimes when the axis of the spring is not aligned with any of the axes of the global coordinate system, two or three pseudo-DOFs per node are assumed for 2D or 3D problems, respectively. In such a case, the two or three DOFs are mathematically related and can be calculated from the angles formed by the axis of the spring and any of the axes of the global coordinate system. Unlike a truss element, in which a constant crosssectional area, a length, and an elastic modulus are required in order for all the entries of the element stiffness matrix to be calculated, only a single spring constant is needed to form the element stiffness matrix of a spring. A cable element can also support an axial load, but it may be further decomposed such that it becomes either a tension-only or compression-only cable element. A tension-only cable element can provide load resistance when it is elongated, but becomes slack and provides no load resistance when it is being shortened. Computationally, this is done by detaching the stiffness when the tension-only cable element goes into compression. A compression-only cable element has an effect opposite to a tension-only cable element. Because a muscle can only generate force when it is in contraction, a compression-only cable element needs to be used to model an active muscle.

2.3 Element

In classical mechanics, a beam is defined as having two of the three dimensions significantly smaller than the third one. The main purpose of a beam is to transmit transverse force through bending. Hence, a beam element has two DOFs per node, a vertical deflection and a rotation, to support both shear force and moment. For this reason, a C1 continuity is required to formulate a classical plane beam element. A beam element is significantly different from a truss element, which supports only axial loading. Another significant difference between a beam and a truss element is that a truss element can either support compression or tension, but not both at the same time. On the other hand, there exists a neutral axis (a plane where all stresses are zero) for a beam element, with one side of the neutral axis subjected to tension while the opposite side is subjected to compression. Although only two nodes are needed to fully describe the geometric orientation of a beam element, some software packages require that a third reference node be added in order to designate the directions for the width and height of the beam. This is useful for proper cross-sectional moments of inertia to be calculated. With this information, the stress within the beam can be calculated at integration points along the direction of the height on the cross section. Other software packages only need the user to directly input the elastic modulus, cross-sectional area, and moment of inertia, and hence no such reference node is needed. In such a case, no stress within the beam would be calculated, due to the lack of cross-sectional shape. A frame element, which can be formulated with the principle of superposition, is in essence the combination of a truss element and a classical plane beam element. This structure is more realistic than classical truss and beam elements, because real-world structures are rarely found to carry purely axial or transverse loading. Many software packages describe a frame element as a general beam element. Modern FEA software packages are versatile and powerful. However, these packages are sometimes so sophisticated, it is difficult to use them properly. In order to create proper models, we need to read the User Manuals as well as the Theoretical Manuals, so the various subtle differences among the structure types can be understood. Table 2.1 summarizes some critical aspects of the aforementioned 1D structure types. Note that these summaries may not be applicable for all FE software packages.

2.3.3 2D ELEMENT TYPE The 2D element type is a 3D structure with one of the dimensions much smaller than the other two (e.g., a plate subjected to in-plane loading). This element type includes, but is not limited to, 3-node linear triangular, 6-node quadratic triangular, 4-node bilinear quadrilateral, 9-node biquadratic quadrilateral, and 8-node serendipity quadrilateral elements (Fig. 2.8). Clearly, the two simplest element types are the 3-node linear triangular and 4-node bilinear quadrilateral elements. The 3-node element type is the easiest one to use for the purpose of meshing.

61

62

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

Table 2.1 One-Dimensional Structure Types That Can Be Modeled With the Simplest 2-Node Element Structure Type

DOFs per Nodea

External Loads Resisted

DOFs

Bar, Truss, Spring Cable

1

Tension and compression

Axial displacement

1

Axial displacement

Beam

2

Frame

5

Either tension or compression Transverse force and bending moment Tension/compression/ bending moment

a

Transverse displacement and rotation Three displacements and two rotations

Pseudo-DOFs used for orientation are not included.

(A)

(D)

(B)

(C)

(E)

FIGURE 2.8 Example 2D element types include: (A) 3-node linear, (B) 6-node quadratic, (C) 4-node bilinear quadrilateral, (D) 9-node biquadratic quadrilateral, and (E) 8-node serendipity quadrilateral elements.

However, many more 3-node elements need to be employed in order to achieve the same accuracy that a 4-node element type can provide. Typically, a constant value is prescribed in the input data deck to describe the thickness (depth) of the entire element. For cases where nonuniform thicknesses are required, some software packages allow different thicknesses to be assigned at

2.3 Element

each node within the element. Importantly, the out-of-plane thickness will not be graphically displayed in any pre- and postprocessing software packages, regardless of the thickness. Allowable limits for the ratio of thickness to in-plane lateral dimensions (width and length) are problem dependent and cannot be definitively prescribed. As a rule of thumb for qualifying the use of this element type, the thickness along the z-axis should be less than 10% of the dimensions in the x- or y-direction. All of these 2D surface element types can be used to represent membrane, plane stress, plane strain, plate, and shell structures. The distinctions among these five structure types are that membrane, plane stress, and plane strain elements possess strength only along the surface of the plane, a plate element has only out-of-plane stiffness (i.e., support load through bending), and a shell element can provide both in-plane and out-of-plane stiffness. All 2D element types can be used independently or in conjunction with 3D element types for modeling 3D structures. For example, a vertebra can be represented by a layer of 2D shell elements to represent the thin outer cortical bone and 3D solid elements to represent the trabecular bone within. The thickness of a membrane element is considered too thin to resist any compressive loading. Despite the thinness, both membrane and plane stress elements can take tensile load. For plane stress elements, any stress components involving an out-of-plane axis must be zero. That is, szz, syz, and szx are all zero for an element that lies on the xey plane. Similarly, a plane strain element in the xey plane postulates that εzz, εyz, and εzx are all zero. In classical mechanics, a shell must have a curved surface (e.g., the outer surface of an oil storage tank or the haul of a submarine), while a plate has a flat surface. However, because of the large number of elements routinely employed in modern FE models, the distinction between curved and flat surfaces has become blurry; a curved surface can be well approximated by joining a large number of flat surfaces. Also, in some software packages, the plate element is allowed to provide stiffness to resist both in-plane and out-of-plane loading. Since the distinction between a shell and a plate is sometimes not as clear in the modern FE method as it was when analytical methods were used, we need to be well aware of the ways plate and shell elements are formulated in our chosen software package, in order to correctly select the proper 2D element type to develop an FE model. Because membrane, plane stress, and plane strain elements provide only in-plane stiffness, all rotational DOFs are constrained, while translational DOFs are allowed. As such, these elements possess two DOFs (two in-plane displacements) per node. Typical applications of membrane elements are related to modeling thin fabrics, such as a seat belt or an airbag. Similarly, the mesothelium, which forms the lining of the pleura (thoracic cavity), pericardium (the membrane enclosing the heart), and peritoneum (abdominal cavity), needs to be modeled using membrane elements. On the other hand, plane stress elements are commonly used for solving problems involving a thin plate subjected to only in-plane loading. For example, a seat-belt buckle, as shown in Fig. 2.9, is commonly modeled using plane stress elements, while a seat belt is commonly modeled with membrane elements.

63

64

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

FIGURE 2.9 A seat-belt buckle modeled using plane stress elements, while a seat belt is modeled with membrane elements.

As explained above, a 2D plate element supports loading that is orthogonal to the plane surface through combined bending and shear. Hence, a 2D plate element has three DOFs per node (a vertical displacement and two in-plane rotations). An example that would be well modeled with plate elements is loading on a cover of a utility access manhole by passing cars. This is because external forces applied to this road structure are, in general, perpendicular to the surface of the cover, which in turn produces bending load. Plate elements can be further categorized into thin and thick plates. Although the distinction between a thin versus a thick plate is not well defined, a thickness to width or length ratio of lower than 10% is generally considered a thin plate, while a ratio greater than 10% is regarded as a thick plate. The Kirchhoff plate theory or KirchhoffeLove plate theory that was derived to calculate the deformations and stresses within a plate is the foundation generally adopted for the formulation of a thin plate. The Love theory (Love, 1888) was developed by Augustus Edward Hough Love (Apr. 1863eJun. 1940), based on a series of lectures given by Gustav R. Kirchhoff (Mar. 1824eOct. 1887) at the University of Heidelberg. Here the birth and death dates are based on the MacTutor History of Mathematics Archive wonderfully compiled by O’Connor and Robertson (2017). Unless otherwise stated, all birth and death dates listed are according to this online source. The assumptions based on Kirchhoff’s proposal require that a straight line that is perpendicular to the midsurface of the plate remains straight and orthogonal to the deformed midsurface (Fig. 2.10). As such, the transverse shear deformation is neglected. For a thick plate, the formulation generally follows the ReissnereMindlin plate theory, which accounts for shear behavior not considered in the Kirchhoff plate theory. The basic differences between Mindlin and Reissner’s plate theories are: •

•

Eric Reissner (Jan. 1913eNov. 1996) assumed that the displacement across the plate (i.e., out-of-plane) may not be linear, and the thickness of the plate may change with loading (Reissner, 1945). Raymond D. Mindlin (Sep. 1906eNov. 1987) assumed a linear variation in the displacement across the plate thickness, and the thickness of the plate remained unchanged with loading. Also, any normal stress through the thickness is neglected (Mindlin, 1951).

2.3 Element

θ

θ

θ=

FIGURE 2.10 A cross-sectional view of a thin plate before and after bending deformation. The Kirchhoff plate theory assumes that a straight line that is perpendicular to the midsurface of the undeformed plate remains straight and perpendicular to the deformed midsurface. The slope or rotation angle is equal to the first derivative of the deflection (i.e., q ¼ dw dx ).

Despite the slight variation between Mindlin and Reissner’s theories, both considered nonzero out-of-plane shear deformations. As such, it is frequently referred to as ReissnereMindlin thick plate theory. This thick plate theory is in contrast to KirchhoffeLove’s thin plate theory, in which zero out-of-plane deformation is assumed. In general, the thick plate formulation is recommended, because this method results in a more accurate solution. However, if there are a large number of elements with high aspect ratios (the ratio of length to width), the thick plate formulation should not be considered. Because the ReissnereMindlin theory is more versatile than other theories, in almost all commercial software packages, such as ABAQUS, ANSYS, LS-DYNA, and PAM-CRASH, the element libraries are based on the ReissnereMindlin theory. In addition to different assumptions used in the derivations of different plate theories, ways of implementing these theories into the element library have also varied. For example, the first type of plate element implemented in LS-DYNA was the Hughes-Liu shell element, while the most computationally efficient type was the Belytschko-Tsay shell element (Hallquist, 2006). As software users, we need to understand which implemented element type best approximates the solution we wish to obtain. A shell element, in essence, is a combination of a plane stress (two DOFs) and a plate element (three DOFs) into one element with five DOFs. Hence, formulating a shell element can be done through superposition of a plane stress element and a plate element. As noted above, in some software packages, deformation of a plate element may include in-plane stretch or shortening and out-of-plane deflection due to bending. To the best of the author’s knowledge, no FE software packages allow users to output deformations due to individual in-plane or out-of-plane loadings. Instead, the in-plane and out-of-plane deformations are summed and reported as the total deformation. Table 2.2 summarizes some critical aspects of the aforementioned 2D structure types.

65

66

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

Table 2.2 Structure Types That Can Be Modeled With the Simplest 4-Node, 2D Surface Element 4-Node Surface Element

Stiffness Matrix

DOF per Node

Available DOF

Membrane, plane stress, and plane strain

2

In-plane translations

88

Plate

3

One out-of-plane translation and two in-plane rotations

1212

Shell

5

Three translations and two in-plane rotations

2020

In the explicit FE solution package LS-DYNA (LSCT, Livermore, CA), the shell element is used for modeling shells, plates, membranes, plane stress, and plane strain, all using the keyword *ELEMENT_SHELL augmented with a large number of element formulation options (ELFORM), such as Belytschko-Tsay membrane (ELFORM 5), plane stress (ELFORM 12), plane strain (ELFORM 13), Belytschko-Wong-Chiang shell (ELFORM 10), etc. As a user, this fact further emphasizes the need to become familiar with the methods of element formulation in software packages before the proper element type is chosen for development of FE models. Due to the limited scope of the current book, additional discussion on this rather complex subject of implementations of various shell element types can

2.4 Formation of Finite Element Mesh

(A)

(B)

(C)

(D)

FIGURE 2.11 Example 3D element types include: (A) 4-node tetrahedral, (B) 8-node trilinear, (C) 6node prismatic, and (D) 5-node pyramidal elements.

not be provided. If you are interested in details on this topic, we recommend reading the LS-DYNA theory manual (Hallquist, 2006).

2.3.4 3D ELEMENT TYPE The three-dimensional (3D) solid element type includes, but is not limited to, 4-node tetrahedral linear, 8-node trilinear, 6-node prismatic or wedge-like, and 5-node pyramidal elements (Fig. 2.11). Higher order elements such as the 10-node quadratic tetrahedral element and the 20-node hexahedron tri-quadratic (serendipity) elements are more complex and not shown in this figure. Each node of the 3D element type possesses only three translational DOFs. Accordingly, moment cannot be applied to this element type, and no rotation will ever occur. Again, the simplest 4-node tetrahedral and 8-node trilinear elements are recommended for the development of 3D solid elements in FE models. The most adaptable element type for any 3D geometric shape is the tetrahedral elements, and these elements are very easy to form using automatic meshing algorithms. However, unless a very large number of elements is used, this element type frequently stiffens the response. On the other hand, high-quality hexahedral elements are more difficult to form, especially for very complex geometry, but usually provide a more accurate solution without a great number of elements.

2.4 FORMATION OF FINITE ELEMENT MESH The two previous sections describe the node and element needed to form the mesh of an FE model. A mesh is a collection of elements interconnected at nodes used to idealize a structure of interest. When a mesh is augmented with proper boundary and loading conditions, it becomes an FE model. For simple problems, such as those highlighted in Fig. 1.14, the mesh can be easily generated by hand, as shown in Table 1.1. A more complex geometry requires the application of an FE preprocessor to develop the mesh. Most modern FE model preprocessors allow users to interchange different coordinate systems (such as the Cartesian, polar, and spherical) for easy development of the model. Nevertheless, the end product is always displayed in Cartesian coordinates. Only a concise definition of the Cartesian coordinate system is listed here. We

67

68

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

recommend that you become familiar with other coordinate systems through personal studies. The other coordinate systems could be useful for identifying the nodal coordinates of objects with special shapes, such as setting nodes around a circle. The Cartesian coordinate system, also known as the rectangular coordinate system, was first used by Rene´ Descartes (Mar. 1596eFeb. 1650). A 2D Cartesian plane consists of two perpendicular axes crossing each other at the origin. In general, the horizontal axis is called the x-axis, and the vertical axis is called the y-axis. For a 3D Cartesian space, a third axis (z-axis) is defined by adding a line through the origin and perpendicular to the xey plane. Using Cartesian coordinates, the x-, y-, and z-coordinates indicate how far away the point is located from the origin. It must be emphasized that the orientation of the element coordinate system is the one that is the most convenient for that element. In other words, the element coordinate system may not need to coincide with the global coordinate system that is used to describe nodal locations. This element-based coordinate system is defined specifically to closely relate to the material behaviors. For example, a femur is a transverse isotropic material that has a higher elastic modulus along its axial direction than on the transverse plane. In this case, a local element coordinate system is needed to define which axis is along the axis of the femur and which other axes define the transverse plane, so that proper material properties can be implemented. Many software packages are available to generate FE meshes. A list of public domain, downloadable, and university developed Automatic Mesh Generation Methods and Software, maintained by a German engineer Robert Schneiders, can be found on the Internet (Schneiders, 2017). Although it is convenient to use these packages to generate FE meshes, we must understand some fundamental background related to the selection of proper element types to correctly develop an FE model. When meshing an FE model, it is essential that any adjacent elements are only connected to neighboring elements through their nodes. In other words, “floating” or “isolated” nodes are not allowed. Fig. 2.12 below demonstrates two neighboring

Floating Node

FIGURE 2.12 Example of a floating node created when two incompatible elements are connected side by side. In this case, continuity for the displacement and strain become questionable.

2.5 Element Shape Functions and [B] Matrix

elements with element 1 on the left side, formed by a 4-node bilinear element, and element 2 on the right side, made of an 8-node serendipity element. We can easily observe that the common edge (the junction between the two elements) consists of two nodes on element 1 and three nodes on element 2. Clearly, there is a floating node. Because element 1 uses a linear interpolation while element 2 uses a quadratic interpolation, the strains may become incompatible along this edge. While some software packages provide warning messages whenever there is a floating node present, other packages may allow users to run the FE model without questioning this potential mistake. Due to the incorrectly formulated input data deck in the latter case, incorrect outputs will be generated.

2.5 ELEMENT SHAPE FUNCTIONS AND [B] MATRIX As briefly mentioned in Section 2.3, the element shape functions are used to interpolate the nodal coordinates to identify the coordinate values of points located anywhere within the element. The same shape functions can also use the modelcalculated nodal results to interpolate physical values located anywhere within the element. Eq. (2.1) explains the concept of the shape functions in mathematical terms, 4x;y;z ¼

n X

Ni 4i

(2.1)

1

where 4 is a physical quantity that could be the nodal coordinate values, nodal temperature, or nodal displacement, 4x,y,z is the 4 value at the coordinate point (x, y, z) within the element, Ni are the shape functions, and 4i are the 4 values at all distinct nodes. The value n is the number of DOFs constituted the element. For example, a 2-node beam element has two DOFs per node for a total of four DOFs, while a 2-node truss element has one DOF per node for a total of two DOFs.

2.5.1 1D, 2-NODE ELEMENT SHAPE FUNCTIONS The simplest 1D element to represent the truss, spring, and cable is the 2-node linear element. A similar 2-node configuration can also be used to represent a beam. However, shape functions for these two categories of elements are very different.

2.5.1.1 2-Node Linear Bar Element Fig. 2.13 shows a 1D element formed by two nodal points, P1 and P2, with displacement at P1 as u1 and displacement at P2 as u2. For a 2-node element, Eq. (2.1) becomes 4x ¼ N1 41 þ N2 42 .

(2.2)

69

70

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

1

FIGURE 2.13 A 2-node linear bar element with nodal displacements u1 at P1 and u2 at P2.

Example 2.1 Assume that the x coordinates for P1 is 4 and P2 is 6. We can intuitively identify the coordinates at one-quarter distance from the left to be 4.5, midpoint to be 5, and three-quarters distance to be 5.5. Obviously, it would be much easier to use a computer program to compute the needed values. The algorithm needed for such a computer program is embedded in the ideas of shape functions. Show the algorithm for finding the shape functions. Solution Assuming a linear interpolation exists within this element, 4 can be expressed as 4x ¼ a0 þ a1 x.

(2.3)

41 ¼ a0 þ a1  4;

(2.4)

42 ¼ a0 þ a1  6:

(2.5)

At P1 (x ¼ 4), and at P2 (x ¼ 6),

To find a1, we subtract Eq. (2.4) from Eq. (2.5). In doing so, we have 4  41 . (2.6) a1 ¼ 2 2 We then insert the results from Eq. (2.6) into either Eq. (2.4) or (2.5). We choose Eq. (2.5): 4  41 4 1 ¼ a0 þ 2  40a0 ¼ 341  242 . (2.7) 2 Finally, we insert a0 and a1 back into Eq. (2.3): 4  41 x. (2.8) 4x ¼ 341  242 þ 2 2 The physical meaning of this equation can be explained as the value of any physical quantity 4 (such as the temperature, nodal coordinate, or nodal displacement) at any location x within the element (between P1 and P2), and it can be calculated using Eq. (2.8). As a trivial example, assuming that the temperatures

2.5 Element Shape Functions and [B] Matrix

are 4 at P1 and 6 at P2, what is the temperature at x ¼ 5? By inserting x ¼ 5, 41 ¼ 4 , and 42 ¼ 6 into Eq. (2.8), we can easily find the result that 4x¼5 ¼ 341  242 þ

42  41 ð6 Þ  ð4 Þ 5 ¼ 5 . x ¼ ð3Þð4 Þ  ð2Þð6 Þ þ 2 2 (2.9)

For a nontrivial example, if the displacements u1 ¼ 41 ¼ 35 mm and u2 ¼ 42 ¼ 47 mm, we can find the displacement at x ¼ 5 from Eq. (2.8) as 4x¼5 ¼ 341  242 þ

42  41 ð47Þ  ð35Þ 5 ¼ 41 mm. x ¼ ð3Þð35Þ  ð2Þð47Þ þ 2 2 (2.10)

The next step seems simple, but it is an essential step is necessary for understanding subsequent material. To find the shape functions, we rearrange Eq. (2.8) by collecting all 41 terms together and all 42 terms together. 4x ¼ ð3  0:5xÞ41 þ ð0:5x  2Þ42

(2.11)

By comparing terms in Eq. (2.11) with those in Eq. (2.2), we can determine that the two shape functions ½ N1 N2  are of the values N1 ¼ 3  0.5x and N2 ¼ 0.5x  2. If we plug the same two 4 values (35 and 47 mm) and x ¼ 5 into Eq. (2.11) 4x¼5 ¼ ð3  0:5  5Þ35 þ ð0:5  5  2Þ47 ¼ 41 mm. We find that the displacement value at the same coordinate is identical to that calculated from Eq. (2.8) and illustrated in Eq. (2.10). The same outcome is easily understandable, because Eq. (2.11) is simply a rearrangement of Eq. (2.8). By writing it in the form of Eqs. (2.2) and (2.11), the element shape functions are easier to visualize. The shape functions shown in Eq. (2.11) are only valid for the 2-node element with P1(x ¼ 4) and P2(x ¼ 6). For different P1 and P2 values, a new set of shape functions are needed. Because deriving shape functions for each element with different coordinate values is time-consuming, the computational cost would be prohibitively high if the FE model were to consist of a great number of elements. A more general approach is to derive the shape functions based on the length of the element. From Fig. 2.13, we define the length of the bar to be L, that is, x(P2)  x(P1) ¼ L.

Example 2.2 Assuming a local coordinate system with x(P1) ¼ 0 and x(P2) ¼ L, use the linear interpolation algorithm for finding the shape functions. Solution Let 4x ¼ a0 þ a1x, as previously shown for the linear interpolation in Eq. (2.3). At P1, (2.12) 41 ¼ a0 þ a1  0; and

71

72

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

at P2, 42 ¼ a0 þ a1  L.

(2.13)

Solving these two simultaneous equations, we have 4  41 . (2.14) a0 ¼ 41 and a1 ¼ 2 L Inserting the values of a0 and a1 into Eq. (2.3), we find that  4  41 x x 4x ¼ 41 þ 2 x ¼ 1  41 þ 42 . (2.15) L L L And the element shape functions ½ N1 N2  are the coefficients listed in Eq. (2.15): N1 ¼

Lx x and N2 ¼ L L

(2.16)

Example 2.3 Similar to the previous two examples, we can assume a different local coordinate L system where xðP1 Þ ¼ L 2 and xðP2 Þ ¼ 2, that is, the new origin of the coordinate system is located at the center of the bar. Determine the algorithm for finding the shape functions. Solution At P1, 41 ¼ a0 þ a1 

L . 2

(2.17)

At P2, L 42 ¼ a0 þ a 1  . 2 Solving these two simultaneous equations, we find that a0 ¼

41 þ 4 2 4  41 and a1 ¼ 2 ; and 2 L

41 þ 42 4 2  4 1 L  2x L þ 2x 41 þ 42 . þ x¼ 2L 2L 2 L Based on Eq. (2.2), the two element shape functions are 4x ¼

N1 ¼

L  2x L þ 2x and N2 ¼ . 2L 2L

(2.18)

(2.19) (2.20)

(2.21)

2.5 Element Shape Functions and [B] Matrix

These three examples show that using different element coordinate systems and different ways to describe the length of the element changes the expressions of the element shape functions. Eq. (2.11) demonstrates that for a bar element with P1 ¼ 4 and P2 ¼ 6, the corresponding shape functions are N1 ¼ 3  0.5x and x N2 ¼ 0.5x  2. From Eq. (2.16), we see that N1 ¼ Lx L and N2 ¼ L for an element with P1 ¼ 0 and P2 ¼ L. Lastly, Eq. (2.21) shows that for a bar element with P1 ¼ L2 and P2 ¼ L2, the corresponding shape functions are N1 ¼ L2x and 2L L þ 2x N2 ¼ 2L . As we can see, element shape functions are different when different coordinate systems are chosen to describe the nodal positions. Now, we observe that the sum of N1 and N2 yields a unit value (1) for all three choices of the corresponding coordinate system. If the point of evaluation is located at P1, we can see that N1 ¼ 1 and N2 ¼ 0 for all three example cases. Similarly, if P2 is selected as the evaluation point, we can see that N1 ¼ 0 and N2 ¼ 1 for all cases. From these exercises, we can deduce the following properties for the element shape functions related to a 2-node bar element: • •

Summation of all the element shape functions is equal to 1. For the two nodal points P1 and P2, N1 ¼ 1 and N2 ¼ 0 when evaluation is at point P1. When the evaluation point is P2, N1 ¼ 0 and N2 ¼ 1.

Later, we can demonstrate that these characteristics also hold true for other element types. An important note is that the sequence of nodal arrangements can greatly affect the outcomes of the calculations for shape functions. In all three examples, the first node is located on the left-hand side of the element, while the second node is located on the right-hand side. Eq. (2.21) is a more general form as compared to Eqs. (2.11) and (2.16), because the origin of its local coordinate system is defined at the center of the element, the same as needed for Gauss quadrature, which is a numerical integration procedure commonly used in the FE method to identify element-related parameters (see Section 4.5.1). Using this coordinate system, the shape functions for different elements of different lengths would still be different. In Section 3.2, the concept of isoparametric shape functions is introduced, which can be used to represent the same type of elements for easy numerical manipulations. Now, we change 4 (a physical quantity) to u (nodal displacement) in Eq. (2.20) and by using the shape functions described in Eq. (2.21), we can write L  2x L þ 2x u1 þ u2 ¼ ½ N1 N2 f u1 u2 gT . (2.22) 2L 2L This equation allows the calculation of the displacement at any point (i.e., u(x)) from nodal displacements u1 and u2, as long as N1 and N2 are known. By definition, the element strain for an axial element can be calculated by differentiating Eq. (2.22) with respect to x: T du df ½ N1 N2 f u1 u2 g ¼ εxx ¼ dx dx

      u1 u1 d½ N1 N2  u1 1 1 ¼ ¼ ¼ ½B . (2.23) dx u2 u2 u2 L L uðxÞ ¼ N1 u1 þ N2 u2 ¼

73

74

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

The difference between nodal displacements f u1 u2 gT and displacements anywhere {u(x)} must be kept in mind to avoid confusion. From Eq. (2.23), we can find the element strain based on the displacements, or more precisely, nodal displacements. Thus, we call this equation the strainedisplacement equation. Along with this equation, we introduce a new [B] matrix (or the strainedisplacement matrix) commonly used in the FE method. This matrix describes the relationship between the strain and nodal displacements and can be directly derived from the shape functions of a bar element by the following equation: d d½N ½ N1 N2  ¼ . (2.24) dx dx In theory, the strainedisplacement [B] matrix should be named strainenodal displacement matrix, while εxx ¼ du dx is better suited for the name straine displacement matrix. The obvious difference arises from the fact that u(x) represents the displacement anywhere within the element, while f u1 u2 gT are nodal displacements. However, the terminology for the [B] matrix has been used universally, and hence we will continue using the term strainedisplacement matrix to represent the [B] matrix. As mentioned, the equation εxx ¼ du dx is, by definition, the straine displacement equation. Because finding the solution to εxx ¼ du dx requires analytical differentiation, it may be difficult to incorporate all needed algorithms into a computer program to do the calculation. Eq. (2.23) implies that the axial strain εxx can be calculated by simple multiplication and summation from the [B] matrix and nodal displacements f u1 u2 gT . Thus, εxx ¼ ½Bf u1 u2 gT is commonly known as the strainedisplacement equation. Therefore, in addition to interpolating nodal values, the element shape functions play another very important role in the FE method, that is, they can be used to calculate the element strain from the [B] matrix and nodal displacements. A final note regarding the constant [B] matrix for a bar element: we can easily deduce from Eq. (2.23) that the axial strain within the 2-node bar element is constant with a 1 magnitude of εxx ¼ u2 u L . ½B ¼

2.5.1.2 2-Node Beam Element Although a beam element uses the same 2-node configuration as the bar, truss, and cable elements, the interpolation functions are quite different. This is because each node of the beam element possesses two DOFs instead of one in any of these elements. The two DOFs are vertical deflection along the z-axis and rotation about the y-axis. Aside from the two nodal points P1 and P2, a third reference node may be needed to prescribe the direction of the vertical z-axis. As mentioned in Section 2.2, this element type requires C1 continuity, and hence shape functions based on Hermite interpolation, which is named after Charles Hermite (Dec. 1822eJan. 1901), are needed. Fig. 2.14 shows a beam element of length L along the x-axis.

2.5 Element Shape Functions and [B] Matrix

θ

θ

FIGURE 2.14 A 2-node beam element of length L along the x-axis has two DOFs (vertical deflection along the z-axis and rotation about the y-axis) per node. The sign convention for a positive moment or rotation is in the counterclockwise direction about the y-axis. Also, a positive force or deflection is along the positive z-axis.

The Hermite interpolation has the form 4ðxÞ ¼ a1 þ a2 x þ a3 x2 þ a4 x3 .

(2.25)

If we let the physical quantity 4(x) represent the vertical deflection w, then wðxÞ ¼ a1 þ a2 x þ a3 x2 þ a4 x3 ¼ N1 w1 þ N2 q1 þ N3 w2 þ N4 q2 .

(2.26)

1

Because of the C continuity, both the vertical deflection (w) and its slope (q, small rotation about the y-axis), which is the first derivative of the deflection, must be continuous. For a very small angle of rotation, qztan q ¼ dw dx , which can be thought of as the slope of w with respect to x. Thus, to find q, we differentiate Eq. (2.26): qðxÞ ¼ a2 þ 2a3 x þ 3a4 x2 .

(2.27)

Applying the boundary conditions yields: wðx ¼ 0Þ ¼ w1 ¼ a1 ;

(2.28)

qðx ¼ 0Þ ¼ q1 ¼ a2 ;

(2.29)

wðx ¼ LÞ ¼ w2 ¼ a1 þ a2 L þ a3 L2 þ a4 L3 ¼ w1 þ q1 L þ a3 L2 þ a4 L3 ; and (2.30) qðx ¼ LÞ ¼ q2 ¼ a2 þ 2a3 L þ 3a4 L2 ¼ q1 þ 2a3 L þ 3a4 L2 .

(2.31)

We write Eqs. (2.30) and (2.31) in matrix form in terms of two unknowns a3 and a4 as " #( ) ( ) a3 w2  w1  q1 L L2 L3 ¼ . 2L 3L2 q2  q1 a4

75

76

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

Using Cramer’s rule, named after Gabriel Cramer (Jul. 1704eJan 1752), the two unknowns a3 and a4 can be easily calculated from Eqs. (2.30) and (2.31) written in matrix form,    w  w  q L L3   2  1 1    q2  q1 3L2  3w2  3w1  2Lq1  Lq2   a3 ¼ ¼ and (2.32)  L2 L3  L2      2L 3L2     L2 w  w  q L   2 1 1     2L  2w1  2w2 þ Lq1 þ Lq2 q2  q1   a4 ¼ ¼ . (2.33)  L2 L3  L3      2L 3L2  For students who are not familiar with Cramer’s rule, consider that we have two simultaneous linear equations: ax þ by ¼ c dx þ ey ¼ f : We can rewrite these two equations in matrix form as

    a b x c ¼ . d e y f Based on Cramer’s rule, the two unknowns f x y gT are solved by finding the ratios of determinants as shown below:     c b a c     f e d f  ; y ¼  ; x¼ a b a b     d e d e where j j represents the determinant of the matrix. More specifically, the denomina

a b tor for both unknowns are the determinant of the matrix . To find the first d e   a unknown x, we make the numerator by replacing the first column of the ded  

a b c nominator matrix with and then calculate the determinant of the d e f new matrix. Similarly, the numerator for the second unknown y is calculated by   c replacing the second column of the matrix with and then computing the f determinant.

2.5 Element Shape Functions and [B] Matrix

Placing the four constants (a1, a2, a3, and a4) back to Eq. (2.26) give us 3w2  3w1  2Lq1  Lq2 2 2w1  2w2 þ Lq1 þ Lq2 3 x þ x . L2 L3 Rearranging the terms in the form of N1w1 þ N2q1 þ N3w2 þ N4q2, as seen in Eq. (2.26), the beam element shape functions are expressed as wðxÞ ¼ w1 þ q1 x þ

wðxÞ ¼

4 X

Ni wi ¼

i¼1

þ

L3  3Lx2 þ 2x3 L3 x  2L2 x2 þ Lx3 3Lx2  2x3 w1 þ q1 þ w2 3 3 L L L3

L2 x2 þ Lx3 q2 . L3 (2.34)

where w(x) is known as the generalized displacement, which includes both the vertical deflection and the rotation. From Eq. (2.34), the four beam element shape functions Ni can be directly visualized. To find the slope q, we need to take the derivative vw ¼ vN1 w þ vN2 q þ vN3 w þ vN4 q . For this reason, we list both N and N for easy i i,x vx 1 vx 1 vx 2 vx 2 vx references. L3  3Lx2 þ 2x3 6Lx þ 6x2 and N ¼ ; 1;x L3 L3

(2.35)

L3 x  2L2 x2 þ Lx3 L3  4L2 x þ 3Lx2 and N2;x ¼ ; 3 L L3

(2.36)

3Lx2  2x3 6Lx  6x2 and N ¼ ; and 3;x L3 L3

(2.37)

N1 ¼ N2 ¼

N3 ¼

L2 x2 þ Lx3 2L2 x þ 3Lx2 and N4;x ¼ . (2.38) 3 L L3 For x ¼ 0 to L, N1 to N4 can be calculated from Eqs. (2.35)e(2.38) and plotted as shown in Fig. 2.15. We can see that all four curves are continuous. We can also envision that the slopes of these four curves are continuous, as well. Hence, this set of shape functions exhibits C1 continuity (i.e., both the deflection and slope of the deflection are continuous). We can see from Fig. 2.15 that N1 ¼ 1, N2 ¼ 0, N3 ¼ 0, and N4 ¼ 0 at x ¼ 0 (i.e., at P1). Similarly, N1 ¼ 0, N2 ¼ 0, N3 ¼ 1, and N4 ¼ 0 can be observed at x ¼ L (i.e., at P2). We can also observe from the plots for N1, N3, and N4 that the slopes (tangential lines) of these three shape-function curves at P1 are zero. Additionally, the slope at P2 is 1 (a positive slope indicates that the rotation is in the counterclockwise direction), as seen from the plot for N2. Similarly, the slope observed from the plot for N4 at P2 is 1, while the slopes seen in the plots at P2 for N1, N2, and N3 are zero. Thus, this set of functions fits the basic characteristics of shape functions. From Eqs. (2.35)e(2.38), we can express the vertical deflection and rotation at any point within the element using the following equations: N4 ¼

wðxÞ ¼ ½ N1 qðxÞ ¼ ½ N1;x

N2 N2;x

N3

N4 f w1

N3;x

q1

N4;x f w1

w2 q1

q2 gT and

(2.39)

T

(2.40)

w2

q2 g .

77

78

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

N1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2 0 0

0 0.2

0.4

0.6

0.8

1

N3

1

0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

-0.2

N4

1 0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2 0 0

N2

1

0 0.2

0.4

0.6

0.8

1

0

0.2

0.4

-0.2

FIGURE 2.15 Shape functions for a 2-node beam element plotted as the ratio of the x coordinates from 0 to L.

Now we shall look into the longitudinal (axial) stress generated during the bending process. Consider an elastic, isotropic, and straight beam subjected to a vertical load at the right end, as shown in Fig. 2.16. From classical mechanics, we know that there exists a neutral axis at which the longitudinal stress vanishes. If the cross section of the beam is symmetric, then the neutral axis is located at the center of the beam. In the configuration shown in Fig. 2.16, any points above the neutral axis (concave side) are in compression while those points below the neutral axis (convex side) are in tension. If we consider the fact that only two nodes are used to represent a beam element, we can intuitively understand that this line element must coincide with the neutral axis of the beam. As such, no longitudinal stress would be observed within this element. Because compressive and tensile stresses described earlier do occur under bending, it would be incorrect not to consider their existence. From Fig. 2.16, the neutral axis is located at half the height of the beam throughout the entire beam. If the bending angle q is small, this angle can be approximated by the slope of the beam, that is, q ¼ dw dx . Now, the axial elongation u at a distance z from the neutral axis in the tensile side can be calculated by multiplying the distance z and bending angle as uðxÞ ¼ z tan qzz  q ¼ z

dw . dx

(2.41)

2.5 Element Shape Functions and [B] Matrix

θ

x

FIGURE 2.16 A beam element with a height of h is bent with a vertical deflection w. It is assumed that a straight line perpendicular to the neutral axis in the undeformed beam remains straight and perpendicular to the deformed neutral axis. The slope of the deformed neutral axis is q ¼ dw dx . Consider a point at a distance z below the neutral axis, the axial elongation (u) at this height is u ¼ z  q, and the maximum axial elongation is located at z ¼ h2 and umax ¼ 12 h  q.

From Eq. (2.41), the longitudinal strainedisplacement equation at a distance z from the neutral axis can now be determined as εxx ¼

du d2 w d2 ¼ z 2 ¼ z 2 ½ N1 dx dx dx

¼ z½Bf w1

q1

w2

N2

N3

N4 f w1

q1

w2

q2 gT

(2.42)

T

q2 g .

The above equation shows the axial strainedisplacement relationship. In other words, the axial strain at a distance z from the neutral axis is the product of z and d 2 w, which is also known as the curvature of the curve. A small curvature indicates dx2 very little change in the bending angle (slope), while a large curvature shows a large, sharp turn. The strainedisplacement (or curvatureedisplacement) matrix [B] for a 2-node beam element can be calculated as shown in the equation below. d2 ½ N1 N2 N3 N4  dx2

6 12x 4 6x 6 12x 2 6x ¼  2þ 3  þ 2 þ   L L L L2 L L L2 L3 ½B14 ¼

(2.43)

2.5.2 2D, 3-NODE LINEAR TRIANGULAR ELEMENT The simplest 2D elements are the 3-node triangular element and 4-node rectangular plane stress element. These elements provide stiffness due to in-plane loading.

2.5.2.1 3-Node Linear Triangular Element Fig. 2.17 shows a triangular element formed by three nodal points, P1, P2, and P3, with two DOFs (translations along the x- and y-axes, i.e., u and v) at each of the nodal points, for a total of six DOFs. Again, the counterclockwise arrangement of the nodal points is followed, as previously prescribed.

79

80

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

,

,

)

) ,

)

FIGURE 2.17 A 3-node triangular element. The coordinates (x,y) and displacements u(x,y) and v(x,y) of any point P within the element can be determined from the respective nodal values.

For the horizontal displacement (i.e., displacement along the x-direction, u(x,y)) at any point (x,y) within the triangle, the following interpolation function can be assumed: uðx; yÞ ¼ a1 þ a2 x þ a3 y ¼ N1 u1 þ N2 u2 þ N3 u3 .

(2.44)

The same relationship can also be assumed for the vertical displacement anywhere within the element: vðx; yÞ ¼ b1 þ b2 x þ b3 y ¼ N1 v1 þ N2 v2 þ N3 v3 .

(2.45)

Based on Eq. (1.10), which is repeated here for convenience 2 3 v 6 07 8 9 6 vx 7 6 7( ) > < εxx > = 6 v7 6 7 u εyy ¼ 6 0 ; 7 vy 7 v > > : ; 6 6 7 gxy 6v v7 4 5 vy vx vv the axial and shear strains can be calculated as εxx ¼ vu vx ¼ a2 ; εyy ¼ vy ¼ b3 ; vu vv and gxy ¼ vy þ vx ¼ a3 þ b2 . As we can see from these calculations, all strain components within the element are constant values. For this reason, the linear 3-node triangular element is also known as the constant strain triangle (CST). Due to this nature, an FE model formed by linear triangular elements typically behaves more stiffly than the physical structure that is being modeled. The only way to avoid this extra stiffness is to provide needed strain variations across the structure by using many elements. For this reason, this type of element is not recommended for complex loading conditions, except for the case where a very fine mesh is used. There is no

2.5 Element Shape Functions and [B] Matrix

universal rule for determining how many triangular elements are needed to achieve an acceptable solution. But it can be seen later from Fig. 2.21 and Table 2.4 in the end of Section 2.5.3 that a simulated beam model made of 1288 CST elements did not perform as well as one made of 40 quadrilateral elements. From Eq. (2.44), three equations for nodal displacements u1, u2, and u3 at the three nodal points P1, P2, and P3, respectively, can be derived as u1 ¼ a1 þ a2 x 1 þ a3 y 1 ;

(2.46)

u2 ¼ a1 þ a2 x2 þ a3 y2 ; and

(2.47)

u3 ¼ a1 þ a2 x3 þ a3 y3 .

(2.48)

Because nodal values u1, u2, and u3 are calculated from FE solutions, the magnitudes of u1, u2, and u3 are considered to be known constant values, not variables. Additionally, x1, x2, and x3 and values y1, y2, and y3 are nodal coordinates, and hence all are known constant values. As such, we have three equations to solve for three unknowns: a1, a2, and a3. There are a number of ways to find these unknowns. Here, the method involving only the matrix operation is discussed. The above three equations can be written in matrix form as 8 9 8 9 2 38a 9 a1 > > > > 1 u 1 x y > > > > > 1 1 < < > = = < 1= 6 7 (2.49) u2 ¼ 4 1 x2 y2 5 a2 ¼ ½F a2 . > > > > > > > > : > ; > : ; ; : u3 1 x3 y3 a3 a3 To find a1, a2, and a3, we can simply invert the [F] matrix and then multiply it to both sides of the equations: 8 9 8 9 > > = < u1 > = < a1 > (2.50) a2 ¼ ½F1 u2 . > > ; : > ; : > a3 u3 As we know from linear algebraerelated courses on matrix operations, the inverse of any matrix [A] can be expressed in terms of its determinant (det) and adjoint matrix (adj) as shown in Eq. (2.51). Note here that a classical adjoint of a square matrix is also known as an adjugate or adjunct of the matrix. ½A1 ¼

1 adjðAÞ detðAÞ

(2.51)

where adjðAÞ ¼ ðcofactor matrix of AÞT .

(2.52)

To find the cofactor matrix of a 33 matrix [A], we assume it has the following form: 2 3 a11 a12 a13 6 7 ½A ¼ 4 a21 a22 a23 5. (2.53) a31 a32 a33

81

82

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

Then, the cofactor matrix of [A] can be expressed as 2 3 A11 A12 A13 6 7 cofactor matrix ðAÞ ¼ 4 A21 A22 A23 5 A31 A32 A33 2       a22 a23   a21 a23   a21      6 þ  þ     6 a   a31   a31 a33 32 a33 6     6   a11 a13   a11 6  a12 a13      6  ¼ 6   þ    a31 a33   a31 6  a32 a33  6      6   a11 a13   a11 6  a12 a13     4 þ    þ  a22 a23   a21 a23   a21

3 a22  7 a32  7 7 7  a12  7 7  7.  a32 7 7 7 a12  7 5 a22 

(2.54)

Note that every other entry has the opposite sign. The process for determining this cofactor matrix is to delete the row and column that is at the position of interest in Eq. (2.54). For example, to find A11, eliminate row 1 and column 1, and then find    a22 a23   with proper sign conventhe determinant of the remaining 22 matrix  a a  32

33

tion. Similarly, we calculate A12 by eliminating row 1 and column 2, finding the determinant of the remaining matrix, and then multiplying by negative 1. Thus,    a21 a23  . Example 2.4 provides an exercise to go through these  A12 ¼  a a  31

33

processes. 2 1 Example 2.4 6 Find the cofactor matrix of A ¼ 4 6 7 Solutions

2 5 8

3

3

7 4 5. 9

      5 4 6 4 6 5       A11 ¼   ¼ 13; A12 ¼   ¼ 26; A13 ¼   ¼ 13 8 9 7 9 7 8       2 3 1 3 1 2       A21 ¼   ¼ 6; A22 ¼   ¼ 12; A23 ¼  ¼6 8 9 7 9 7 8       2 3 1 3 1 2       A31 ¼   ¼ 7; A32 ¼   ¼ 14; A33 ¼   ¼ 7 5 4 6 4 6 5 2 3 13 26 13 6 7 Thus, cofactor matrixðAÞ ¼ 4 6 12 6 5. 7 14 7

2.5 Element Shape Functions and [B] Matrix

Next, we find the adjoint by transposing the cofactor matrix based on Eq. (2.52). 3 2 A11 A21 A31 7 6 7 adjðAÞ ¼ ðcofactor ðAÞÞT ¼ 6 4 A12 A22 A31 5 A13 A23 A31 2       3 2     a22 a23   a13 a12   a22 a23   a12 a13   a12 a13            6       7 6    6 6 a 7         a32 a33 a22 a23 7 6 a32 a33 a33 a32  6 32 a33 6      7 6    6 6 a  a11 a13  7 6  21 a23   a11 a13  7 6  a23 a21   a11 a13    ¼6 ¼6     7    6 6  a 7 6 a31 a33   a31 a33  7 6  a33 a31   a31 a33  a 21 23 6 7 6      7 6    6   a11 a12   a11 a12  7 6  a21 a22   a12 a11  6  a21 a22            4 5 4           a         a a a a a a a a a  31

32

31

32

21

22

31

32

32

31

(2.55) The last part of Eq. (2.55) is aimed at eliminating the minus signs by swapping columns. However, there is no need to do so, except for the cosmetic reasons. To explain this procedure, consider that the determinant for the matrix 

  a b     c d  ¼ ad  bc. By swapping the first and second columns,   b a    d c  ¼ ðad  bcÞ. Thus, swapping columns results in negative determinants, which in turn allows the removal of the negative signs. From Eq. (2.50), the three constants a1, a2, and a3 can be found by using the following equation: 8 9 8 9 2 31 8 9 1 x1 y1 > > > < a1 > = < u1 > = < u1 > = 6 7 (2.56) a2 ¼ ½F1 u2 ¼ 4 1 x2 y2 5 u2 . > > > : > ; : > ; : > ; a3 u3 1 x3 y3 u3 From Eqs. (2.51) and (2.55), we find 2  F22 6  6F 6 32 6 6  F23 1 1 6 ½F1 ¼ adjðFÞ ¼ 6 detðFÞ detðFÞ 6  F33 6 6 6  F21 4  F31

 F23   F33   F21   F31   F22   F32 

  F13    F33   F11    F31   F12    F32

 F12   F32   F13   F33   F11   F31 

  F12    F22   F13    F23   F11    F21

3 F13  7 F23  7 7 7  F11  7 7  7. F21  7 7 7 F12  7 5 F22  (2.57)

83

  a12   a 22   a13   a 23   a11   a 21

3 a13  7 a23  7 7 7 a11  7 7 7 7 a21  7 7 7 a12  7 5  a  22

84

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

Replacing the known F11 to F33 values in Eq. (2.57) results in 2     3  x2 y2   y1 x1   x1 y1       7 6 6  x y   y x   x y  7 3 3 2 2 7 6 3 3     7 6      6 1 6  y2 1   1 y1   y1 1  7 7 1 ½F ¼      7. 6 detðFÞ 6  y3 1   1 y3   y2 1  7 6   7   7 6 6  1 x2   x1 1   1 x1  7 4   5    1 x3   x3 1   1 x2 

(2.58)

Finally, we simplify the matrix by solving each of the determinants as follows: 2 3 x2 y3  x3 y2 x3 y1  x1 y3 x1 y2  x2 y1 1 6 7 ½F1 ¼ (2.59) y3  y1 y1  y2 5 4 y2  y3 detðFÞ x3  x2 x1  x3 x2  x1 We know that the determinant of the [F] matrix shown in Eq. (2.59) is related to the area of the triangle. The processes involved in deriving the formula are relatively lengthy and will not be covered in this book. However, the results are as follows: 2 3 1 x1 y1 1 1 6 7 Area of a triangle ¼ A ¼ det½F ¼ det4 1 x2 y2 5. (2.60) 2 2 1 x3 y3 From Eqs. (2.56), (2.59) and (2.60), we have 8 9 8 9 u1 > a1 > > > > > > > > > > > > > = < > = < > 1 a2 ¼ ½F u2 > > > > > > > > > > > > > > ; : > ; : > a3 u3 2 x2 y3  x3 y2 x3 y1  x1 y3 6 6 1 6 y2  y3 ¼ y3  y1 detðFÞ 6 4 x3  x2 x1  x3 2 ¼

1 6 4 2A

x2 y3  x3 y2

x3 y1  x1 y3

y2  y3

y3  y1

x3  x2

x1  x3

ðEq. 2:56Þ

38 9 x1 y2  x2 y1 > > u1 > > > > > 7> 7< = 7 ðEq. 2:59Þ y1  y2 7 u2 > > > 5> > > > : > ; x2  x1 u3

38 9 x1 y2  x2 y1 > = < u1 > 7 y 1  y 2 5 u2 > ; : > x2  x1 u3

ðEq. 2:60Þ

2.5 Element Shape Functions and [B] Matrix

Inserting a1, a2, and a3 into Eq. (2.44) and then rearranging terms in the form of u(x) ¼ N1u1 þ N2u2 þ N3u3, we have uðxÞ ¼ a1 þ a2 x þ a3 y ¼

¼

1 ½ðx2 y3  x3 y2 Þu1 þ ðx3 y1  x1 y3 Þu2 þ ðx1 y2  x2 y1 Þu3  2A

þ

1 ½ðy2  y3 Þu1 þ ðy3  y1 Þu2 þ ðy1  y2 Þu3 x 2A

þ

1 ½ðx3  x2 Þu1 þ ðx1  x3 Þu2 þ ðx2  x1 Þu3 y 2A

1 ½ðx2 y3  x3 y2 Þ þ ðy2  y3 Þx þ ðx3  x2 Þyu1 2A þ

1 ½ðx3 y1  x1 y3 Þ þ ðy3  y1 Þx þ ðx1  x3 Þyu2 2A

þ

1 ½ðx1 y2  x2 y1 Þ þ ðy1  y2 Þx þ ðx2  x1 Þyu3 2A

Compared with Eq. (2.44), the three shape functions of a linear triangular element can be expressed as N1 ¼

1 f ðx2 y3  x3 y2 Þ þ ðy2  y3 Þx þ ðx3  x2 Þyg; 2A

(2.61)

1 (2.62) f ðx3 y1  x1 y3 Þ þ ðy3  y1 Þx þ ðx1  x3 Þyg; and 2A 1 (2.63) N3 ¼ f ðx1 y2  x2 y1 Þ þ ðy1  y2 Þx þ ðx2  x1 Þyg. 2A As previously mentioned, this same set of shape functions can be used to find the coordinates, horizontal displacement u(x,y), and vertical displacement v(x,y) of any point P(x,y) within the element. The three shape functions have geometric meanings. Fig. 2.17 shows a point P(x,y) that divides the whole triangle into three subtriangles with areas A1, A2, and A3 opposite to points P1, P2, and P3, respectively. From Eq. (2.60), 3 2 1 x y 7 6 2A1 ¼ det4 1 x2 y2 5 ¼ ðx2 y3  x3 y2 Þ þ ðy2  y3 Þx þ ðx3  x2 Þy. (2.64) N2 ¼

1

x3

y3

Note that Eqs. (2.64) and (2.61) have the same attributes. In other words, A1 . A Similarly, N2 and N3 can be expressed in terms of A2 and A3 as N1 ¼

N2 ¼

A2 A3 and N3 ¼ . A A

(2.65)

(2.66)

85

86

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

FIGURE 2.18 Three area coordinates N1, N2, and N3 of a triangular element. We can see from this figure that N1 þ N2 þ N3 ¼ 1 holds true at any point within the element.

These three equations demonstrate that the shape functions for a 3-node linear constant strain triangle are related to the ratios of the three subareas to the entire area of the triangle. Hence, these shape functions are referred to as the area coordinates. Note that there are only three shape functions, and together they comprise the entire area of the triangle. Fig. 2.18 graphically shows the magnitudes of these three area coordinates at different locations within the triangle. The strainedisplacement equations {ε} ¼ [B]{u} for this triangular element can be derived from u and v using Eq. (1.10) restated for a triangular element as listed below.

εxx ¼

uðx; yÞ ¼ N1 u1 þ N2 u2 þ N3 u3

(2.67)

vðx; yÞ ¼ N1 v1 þ N2 v2 þ N3 v3

(2.68)

vu vN1 vN2 vN3 ¼ u1 þ u2 þ u3 vx vx vx vx

1 ½ ðy2  y3 Þu1 þ ðy3  y1 Þu2 þ ðy1  y2 Þu3  ¼ 2A εyy ¼

vv vN1 vN2 vN3 ¼ v1 þ v2 þ v3 vy vy vy vy

1 ½ ðx3  x2 Þv1 þ ðx1  x3 Þv2 þ ðx2  x1 Þv3  ¼ 2A gxy ¼ ¼

(2.69)

(2.70)

vu vv vN1 vN2 vN3 vN1 vN2 vN3 þ ¼ u1 þ u2 þ u3 þ v1 þ v2 þ v3 vy vx vy vy vy vx vx vx 1 ½ ðx3  x2 Þu1 þ ðy2  y3 Þv1 þ ðx1  x3 Þu2 þ ðy3  y1 Þv2 2A

þ ðx2  x1 Þu3 þ ðy1  y2 Þv3 

(2.71)

2.5 Element Shape Functions and [B] Matrix

In matrix form, the strainedisplacement equations and the strainedisplacement matrix [B] are expressed as 8 9 2 3 ðy2  y3 Þ 0 ðy3  y1 Þ 0 ðy1  y2 Þ 0 > < εxx > = 1 6 7 εyy ¼ 0 ðx3  x2 Þ 0 ðx1  x3 Þ 0 ðx2  x1 Þ 5 4 > > 2A : ; gxy ðx3  x2 Þ ðy2  y3 Þ ðx1  x3 Þ ðy3  y1 Þ ðx2  x1 Þ ðy1  y2 Þ 8 9 u1 > > > > > > > > > > > > v 1 > > > > > > > = 2 .  > v2 > > > > > > > > > > > u3 > > > > > > > > : ; v3

(2.72)

Recalling that the strainedisplacement matrix [B] describes the relationship between the8 strain tensor and nodal displacements, we write 9 u > > 1 > > > > > > > > v1 > > 8 9 > > > > > > > > < εxx > < u2 > = = εyy ¼ ½B , and therefore v2 > > > > : > > ; > > gxy > > u3 > > > > > > > > > > > > : ; v3 2 3 ðy2  y3 Þ 0 ðy3  y1 Þ 0 ðy1  y2 Þ 0 1 6 7 ½B ¼ 0 ðx3  x2 Þ 0 ðx1  x3 Þ 0 ðx2  x1 Þ 5. 4 2A ðx3  x2 Þ ðy2  y3 Þ ðx1  x3 Þ ðy3  y1 Þ ðx2  x1 Þ ðy1  y2 Þ (2.73)

2.5.3 4-NODE RECTANGULAR BILINEAR PLANE ELEMENT WITH EDGES PARALLEL TO THE COORDINATE AXES Fig. 2.19 shows a 4-node rectangular plane element formed by points P1, P2, P3, and P4, with the origin of the coordinate system located at the geometric center of the element. Again, these four nodes are arranged in a counterclockwise manner, as described earlier. Each node has two DOFs, translations along the x-axis (u) and along the y-axis (v). Further, it is assumed that the same shape functions are used to interpolate nodal coordinates, horizontal displacement (u), and vertical displacement (v), at any point within the element from corresponding nodal values. We assume that the following polynomials are used to interpolate u and v: uðx; yÞ ¼ a1 þ a2 x þ a3 y þ a4 xy ¼ N1 u1 þ N2 u2 þ N3 u3 þ N4 u4 and

(2.74)

vðx; yÞ ¼ b1 þ b2 x þ b3 y þ b4 xy ¼ N1 v1 þ N2 v2 þ N3 v3 þ N4 v4 .

(2.75)

87

88

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

FIGURE 2.19 A 4-node rectangular bilinear element with dimensions 2a  2b with the origin of the coordinate system located at the center of the element.

where a1, a2, a3, and a4 and b1, b2, b3, and b4 are constants, and N1, N2, N3, and N4 are the element shape functions. Note that the interpolation equations involving constants a’s and b’s have the same form regardless of which origin of the Cartesian coordinate system chosen. However, selecting a different coordinate system would result in a different set of shape functions, as noted in the descriptions for bar elements (see Section 2.5.1). This element type is called the bilinear element because the u and v displacements are linear along the x-axis when a constant y-coordinate is chosen and linear along the y-axis when a constant value of the x-coordinate is selected. For nodal displacements u1 through u4, the following equations can be written from Eq. (2.74) by plugging in the respective nodal coordinates: u1 ¼ uða; bÞ ¼ a1 þ a2 ðaÞ þ a3 ðbÞ þ a4 ab

(2.76)

u2 ¼ uða; bÞ ¼ a1 þ a2 ðaÞ þ a3 ðbÞ þ a4 ðabÞ

(2.77)

u3 ¼ uða; bÞ ¼ a1 þ a2 ðaÞ þ a3 ðbÞ þ a4 ab

(2.78)

u4 ¼ uða; bÞ ¼ a1 þ a2 ðaÞ þ a3 ðbÞ þ a4 ðabÞ

(2.79)

The four unknowns (a1, a2, a3, and a4) can be found by solving these four 3 þu4 2 þu3 u4 simultaneous equations. These values, a1 ¼ u1 þu2 þu , a2 ¼ u1 þu4a , 4  u1  u2 þ u3 þ u4 u1  u2 þ u3  u4 , and a4 ¼ , are then inserted back into Eq. (2.74). a3 ¼ 4b 4ab Then we rearrange all terms in the order of u1, u2, u3, and u4, so that horizontal displacement at any point, u(x,y), within the element can be found using the following equation: uðx; yÞ ¼ N1 u1 þ N2 u2 þ N3 u3 þ N4 u4 ¼

ða  xÞðb  yÞ ða þ xÞðb  yÞ u1 þ u2 4ab 4ab þ

ða þ xÞðb þ yÞ ða  xÞðb þ yÞ u3 þ u4 . 4ab 4ab

(2.80)

2.5 Element Shape Functions and [B] Matrix

That is, the shape functions for a 4-node bilinear element are displayed as ða  xÞðb  yÞ 4ab ða þ xÞðb  yÞ N2 ¼ 4ab (2.81) ða þ xÞðb þ yÞ N3 ¼ 4ab ða  xÞðb þ yÞ . N4 ¼ 4ab Because Eqs. (2.74) and (2.75) are of the same form, we can easily deduce that the same set of shape functions can be used to find the vertical displacement v(x,y) at any of the points within the element. Again, it can be checked that this set of shape functions fits the criteria that: (1) the sum of all shape functions is equal to 1 and (2) at point P1(a,b), N1 ¼ 1, while N2 N3 ¼ N4 ¼ 0, etc. We can now derive the strainedisplacement [B] matrix from Eqs. (2.74), (2.75), and (2.81) as N1 ¼

3 2 3 8P 9 4 v v 6 0 7> > > Ni ui > 07 > > 6 7 6 vx > 8 9 6 vx > > 7 7> 6 i¼1 < = ε 8 9 > 6 7 7 xx > 6 > > > > 7 > > 4 7 7>X 6 εyy ¼ 6 > > > 6 0 vy 7 : ; ¼ 6 0 vy 7 > > > > > N v : > > 6 7 v 7 6 i i; > > : ; 6 7 7 i¼1 6 gxy 6 7 7 6 4v v5 6v v7 5 4 vy vx vy vx 2

2

vN 6 1 6 vx 6 6 6 ¼6 0 6 6 6 vN1 4 vy

2 1 6 ¼ 4 4ab

0

vN2 vx

0

vN3 vx

0

vN4 vx

vN1 vy

0

vN2 vy

0

vN3 vy

0

vN1 vx

vN2 vy

vN2 vx

vN3 vy

vN3 vx

vN4 vy

ðb  yÞ

0

ðb  yÞ

0

ða  xÞ

0

ða  xÞ ðb  yÞ ða þ xÞ

8 9 u1 > > > > > > > > > > > > v 1 > > 3> > > > > > > > > > > > u2 > > 0 7> > > 7> > > 7 < v2 = vN4 7 7> > 7 > u3 > vy 7 > > > > > 7> > > > > > v3 > > > vN4 7 5> > > > > > vx > > > u4 > > > > > > : > ; v4

8 9 u1 > > > > > > > > > > > v1 > > > > > > > > > > > > > u > > 2 > > > > > > 3> > < v2 = 0 ðb þ yÞ 0 ðb þ yÞ 0 7> > ða þ xÞ 0 ða þ xÞ 0 ða  xÞ 5 > u > > > 3> > > > > ðb  yÞ ða þ xÞ ðb þ yÞ ða  xÞ ðb þ yÞ > > > > > v3 > > > > > > > > > > u4 > > > > > > > : > ; v4

89

90

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

2 1 6 ½B ¼ 4 4ab

ðb  yÞ

0

ðb  yÞ

0

ðb þ yÞ

0

ðb þ yÞ

0

ða  xÞ

0

ða þ xÞ

0

ða þ xÞ

0

ða  xÞ ðb  yÞ ða þ xÞ

ðb  yÞ

ða þ xÞ ðb þ yÞ

ða  xÞ

0

3

7 ða  xÞ 5. ðb þ yÞ (2.82)

Another commonly used method to obtain the element shape functions is the application of the Lagrange interpolation, published by Joseph-Louis Lagrange (Jan. 1736eApr. 1813). According to Wolfram Math World (2017), the Lagrange interpolation was first published by Waring in 1779, rediscovered by Euler in 1783, and then published by Lagrange in 1795 (Jeffreys and Jeffreys, 1988). Using this interpolation method, a polynomial curve of up to (n  1)th degrees can be found to pass through a set of n data points on the xey plane. For example, a set of three data points will need a second degree polynomial to pass through all three points. Assume a set of n data points are located at (x1, y1), (x2, y2), ., (xn, yn) on the xey plane, described as yi ¼ f(xi). Then, the Lagrange interpolation of this function has the form of f ðxÞ ¼ Q1 y1 þ Q2 y2 þ . þ Qn yn ; where

ðx  x2 Þðx  x3 Þ.ðx  xn Þ ; ðx1  x2 Þðx1  x3 Þ.ðx1  xn Þ

(2.84)

ðx  x1 Þðx  x3 Þ.ðx  xn Þ ; and ðx2  x1 Þðx2  x3 Þ.ðx2  xn Þ

(2.85)

ðx  x1 Þðx  x2 Þ.ðx  xn1 Þ . ðxn  x1 Þðxn  x2 Þ.ðxn  xn1 Þ

(2.86)

Q1 ¼ Q2 ¼

(2.83)

Qn ¼

Note that if x ¼ x1, the numerator and denominator are the same for Eq. (2.84), and thus Q1 ¼ 1. All other equations have x  x1 in the numerator, which makes Q2 through Qn equal to 0. The same findings also apply to x ¼ x2, x ¼ x3, ., x ¼ xn. Finally, the sum of Q1 through Qn equals 1. Thus, functions Q1 through Qn fit the characteristics of the shape functions. For this reason, the Lagrange interpolation method is frequently used to define the element shape functions, especially for those high-order elements. Example 2.5 Find a polynomial equation f(x) that passes through three points on the xey plane with coordinates of (1, 1), (2, 4), and (3, 9). Solution Because only three points are involved, Eqs. (2.84)e(2.86) become ðx  x2 Þðx  x3 Þ Q1 ¼ ðx1  x2 Þðx1  x3 Þ Q2 ¼

ðx  x1 Þðx  x3 Þ ðx2  x1 Þðx2  x3 Þ

Q3 ¼

ðx  x1 Þðx  x2 Þ ðx3  x1 Þðx3  x2 Þ

2.5 Element Shape Functions and [B] Matrix

We plug the three points into Eqs. (2.84)e(2.86): Q1 ¼

ðx  2Þðx  3Þ x2  5x þ 6 ¼ ; ð1  2Þð1  3Þ 2

Q2 ¼

ðx  1Þðx  3Þ x2  4x þ 3 ¼ ; and ð2  1Þð2  3Þ 1

Q3 ¼

ðx  1Þðx  2Þ x2  3x þ 2 ¼ . ð3  1Þð3  2Þ 2

We now have our three Q values needed for Eq. (2.83), and the y values are y1 ¼ 1, y2 ¼ 4, and y3 ¼ 9. By plugging these three values into Eq. (2.83), we find that x2  5x þ 6 x2  4x þ 3 x2  3x þ 2 1þ 4þ 9 2 1 2

  1 9 5 27 ¼  4 þ x2 þ  þ 16  x þ ð3  12 þ 9Þ ¼ x2 . 2 2 2 2

f ðxÞ ¼

For each x value on a polynomial equation, there is only one corresponding y value. As such, x1 s x2 s x3.sxn as long as there are no repeated points on the polynomial equation. Hence, it is not possible for the denominator to become zero, and there is no need to worry about the problem associated with division to zero.

To find the shape functions for a 2D, 4-node plane element, we need to employ the Lagrange interpolation twice, first along the x-direction and then along the y-direction, or vice versa. Finally, we integrate the shape functions obtained from each axis into one set of shape functions for this 2D element. The following section shows a step-by-step approach to determine the 2D, 4-node plane element shape functions. Step 1: We first determine the two shape functions to be used for interpolating any physical value 4 from left to right. Considering x1 ¼ a and x2 ¼ a, the two shape functions needed for interpolation along the x-direction are N1x ¼

x  x2 xa xa x  x1 xþa xþa ¼ and N2x ¼ ¼ ; ¼ ¼ 2a 2a x1  x2 a  a x2  x1 a þ a

(2.87)

where the “1” in N1x represents node 1 and the “2” in N2x represents node 2. Step 2: Next, we determine the two shape functions to be used for interpolating any physical value 4 from bottom to top. Considering y1 ¼ b and y4 ¼ b, the two shape functions needed for interpolation along the y-direction are N1y ¼

y  y2 yb yb y  y1 yþb yþb ¼ and N4y ¼ ¼ ; ¼ ¼ 2b 2b y1  y2 b  b y4  y1 b þ b

(2.88)

where the “1” in N1y represents node 1, and the “4” in N4x represents node 4.

91

92

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

Step 3: For a line parallel to the x-axis at y ¼ b (i.e., line P1eP2 in Fig. 2.19), a physical value anywhere on this line (4a) can be written using the shape functions N1x and N2x and the physical values at node 1 (41) and node 2 (42) as xa xþa 41 þ 4 . (2.89) 2a 2a 2 We can use the same two shape functions N1x and N2x to interpolate the physical values for any points on the line P4eP3 (4b) using the physical values at node 4 (44) and node 3 (43) as 4a ¼ 4at any point on P1 P2 ¼ N1x 41 þ N2x 42 ¼

4b ¼ 4at any point on P4 P3 ¼ N1x 44 þ N2x 43 ¼

xa xþa 44 þ 4 . 2a 2a 3

(2.90)

Step 4: Next, we interpolate from 4a (corresponding to the bottom line of the rectangle, i.e., line P1eP2) to 4b (corresponding to the top line of the rectangle, i.e., line P4eP3) along the y-axis so that the physical value for any point within the element can be obtained:

  yb xa xþa yþb xa xþa 41 þ 42 þ 44 þ 43 4 ¼ N1y 4a þ N4y 4b ¼ 2b 2a 2a 2b 2a 2a ¼

ðy  bÞðx  aÞ ðy  bÞðx þ aÞ ðy þ bÞðx  aÞ ðy þ bÞðx þ aÞ 41  42  44 þ 43 . 4ab 4ab 4ab 4ab (2.91)

Step 5: Finally, we organize the terms in Eq. (2.91), so the x is the second component in the first set of parentheses, and y is the second component in the second set. Additionally, we rectify the negative signs associated with 42 and 44. After these procedures, the resulting four shape functions are N1 ¼

ða  xÞðb  yÞ ; 4ab

(2.92)

N2 ¼

ða þ xÞðb  yÞ ; 4ab

(2.93)

ða þ xÞðb þ yÞ ; and 4ab

(2.94)

N3 ¼

ða  xÞðb þ yÞ . (2.95) 4ab This set of four shape functions is identical to the one shown in Eq. (2.81), which was found using a different approach. In addition to the aforementioned applications, shape functions are useful for calculating distributions of any physical quantities, such as the temperature contours on a weather map. N4 ¼

2.5 Element Shape Functions and [B] Matrix

Example 2.6 The weather report shows that the temperatures in the cities of Novi, Royal Oak, Troy, and Commerce Township are 60, 64, 68, and 62 F, respectively. Assume that the coordinates, in miles, of these four cities are (10, 5), (10, 5), (10, 5), and (10, 5), respectively. Draft the temperature contours based on the locally measured temperatures and geometric locations based on bilinear interpolation. Solution Based on the coordinates of the four cities, it is clear that these four cities are located at the four corners of a rectangle, and the origin of the coordinate system is located exactly at the center of the rectangle. Thus, the rectangle formed by the four cities is analogous to the rectangle shown in Fig. 2.19, with a ¼ 10 and b ¼ 5. From Eqs. (2.92)e(2.95), we find that ð10  xÞð5  yÞ ð10 þ xÞð5  yÞ ð10 þ xÞð5 þ yÞ ; N2 ¼ ; N3 ¼ ; N4 200 200 200 ð10  xÞð5 þ yÞ . ¼ 200

N1 ¼

The physical values (temperatures) at the four corners are 41 ¼ 60 F, 42 ¼ 64 F, 43 ¼ 68 F, and 44 ¼ 62 F and the interpolation equation is T(x,y) ¼ N141 þ N242 þ N343 þ N444. Combining the values with the equation gives us ð10  xÞð5  yÞ ð10 þ xÞð5  yÞ ð10 þ xÞð5 þ yÞ  60 þ  64 þ  68 Tðx; yÞ ¼ 200 200 200 þ

ð10  xÞð5 þ yÞ  62 200

ð50  5x  10y þ xyÞ  60 þ ð50 þ 5x  10y  xyÞ  64 þð50 þ 5x þ 10y þ xyÞ  68 þ ð50  5x þ 10y  xyÞ  62 ¼ 200 ¼

12; 700 þ 50x þ 60y þ 2xy 200

To draft the contours, let us first identify locations where the temperature is the same based on the above equation. We start with 62 F. 12; 700 þ 50x þ 60y þ 2xy ¼ 62 050x þ 60y þ 2xy ¼ 300 200 Because we have only one equation to solve for two unknowns, multiple combinations of (x, y) could satisfy the equation listed above. Rewrite this equation in the form of ð60 þ 2xÞy ¼ 300  50x0y ¼

300  50x . 60 þ 2x

93

94

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

Using this equation, we can assume an x value within the ranges of 10  x  10, and then calculate the corresponding y value. If the calculated y value is outside the boundary of 5  y  5, then this point is not valid.

FIGURE 2.20 The left figure shows several northwestern suburbs of Detroit based on a Google map. Assume that the four cities are at the four corners of a rectangle of 20  10 miles in size. The 4-node plane shape functions can be used to find the isotherm lines shown on the right, based on locally measured temperatures at the corner cities.

Table 2.3 Geometric Points for 62 F Isotherm x

10.0

9.0

8.0

7.0

6.0

5.0

4.0

3.0

2.0

1.0

0.0

y

5.0

3.57

2.27

1.09

0.0

1.0

1.92

2.78

3.57

4.31

5.0

For illustration purposes, we assume that x ¼ 10, 9, ., 0, 1, ., 10 and then calculate the corresponding y values. For an isotherm of 62 F, locations listed in Table 2.3 satisfy the equation above. Also, for any locations with x > 0, the corresponding y values are out of the range and cannot be used. We repeat the procedure to find the geometric points associated with the 64 and 66 F isotherm lines. Connecting these points for the isotherms yields the contour map shown on the right of Fig. 2.20.

2.5.3.1 Comparison of CST and Bilinear Quadrilateral Element

Example 2.7 Assume a cantilever beam of length 1 m is loaded at the free end by a downward vertical force of 1000 N (Fig. 2.21 top). The beam has a Young’s modulus of 100 GPa, a cross-sectional width of 0.05 m, and height of 0.1 m, which result

2.5 Element Shape Functions and [B] Matrix

FIGURE 2.21 A cantilever beam (top) is fixed on the left side and vertically loaded on the right side. Two example FE models consisting only of triangular and quadrilateral plane elements are presented to represent the beam. The model made of triangular elements has a total of 82 elements (middle), while the model made of quadrilateral elements has a total of 640 elements (bottom). 5  105 m4 . In theory, a beam should be in a cross-sectional moment of inertia of 12 modeled with beam elements to secure the best result. However, in certain circumstances plane elements integrate better with other elements within the same structure. Calculate the free-end deflection (maximum deflection) using (1) an analytical approach and (2) FE models made of 2-node beam elements, 3-node constant strain triangular elements, and 4-node bilinear quadrilateral elements of varying element sizes. Calculate the maximum deflection of each model using an FE software package and tabulate the results.

Solution For a cantilever beam of length L, elastic modulus E, moment of inertia I, and load P at the free end (Fig. 2.21 top), the maximum deflection based on analytical solution is PL3/3EI. Thus, the exact solution results in the maximum deflection of 0.8 mm. Two example cantilever beam FE models consisting of a total of 82 2D, triangular surface elements and 640 2D, quadrilateral surface elements are shown in the middle and bottom of Fig. 2.21, respectively. Additional models based on 2-node beam elements, 3-node CST elements, and 4-node quadrilateral elements have also been created, as listed in Table 2.4. An FE software package ANSYS 14.5 (Canonsburg, PA) was used to calculate the maximum deflections for all

95

96

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

Table 2.4 Comparisons of Maximum Deflections Calculated by an FE Model Against Analytical Analysis Element Type 2-Node beam element

3-Node CST plane element 4-Node bilinear plane element

No. of Row/Element 1/2 1/4 1/8 2/82 4/324 8/1288 2/40 4/160 8/640

Max. Deflection (mm) 0.765823 0.796087 0.803653 0.537724 0.716967 0.782914 0.805085 0.806810 0.807759

Difference (%) 4.27 0.49 0.46 32.78 10.34 2.14 0.64 0.85 0.97

Models consisting of beam elements are very efficient, as demonstrated by the fact that only four elements are needed to achieve an accuracy within 0.5%, while models made of triangular elements require many more elements to match accuracy obtained when other element types are selected.

models, and results are compared to the analytical solutions (Table 2.4). From these results, it is clear that the beam element is the most efficient one to use. In theory, a structure consists of only beam elements and modeling with only beam elements should provide the same result as that calculated using analytical method. However, roundoff error and the number of integration points selected contribute to the small error observed. Nevertheless, the beam element is the best element type to choose for modeling beam structures. Also, it can be seen that the model based on beam elements converges much quicker than the model that consists of CST or bilinear elements. For models meshed entirely with CST elements, the error is greater than 2%, even with a fine mesh of 1288 elements. This element type converges slowly, and hence a model of this type requires a large number of elements to predict acceptable solutions. On the other hand, quadrilateral elements converge much quicker, and fewer elements are needed to attain the same results as what can be achieved using the CST element type. This simple example demonstrates that the CST element should be used sparsely.

2.5.4 2D, 4-NODE PLATE ELEMENT SHAPE FUNCTIONS WITH EDGES PARALLEL TO THE COORDINATE AXES 2.5.4.1 Use Pascal’s Triangle to Select Polynomial Terms A plate is actually a 3D structure component that is typically simplified using the 2D element type. As mentioned in Section 2.3, a plate element possesses three DOFs, a

2.5 Element Shape Functions and [B] Matrix

vertical deflection (w) and two in-plane rotations (qx and qy), per node for a total of 12 DOFs. Instead of having a separate plate element type, some software packages simply use a “generalized shell element” with several options for representing the membrane element, plate element, and shell element. Considering the lumping of these three element types together into one generalized shell element, there can be a great deal of confusion for software users. In addition, significant differences exist in the shell element formulations among different software packages. It is therefore necessary for users of software to read and comprehend the contents discussed in the theory and user manuals published by the software vendors. In other words, it is in the user’s best interest to understand the fundamental theory behind the formulation of each element type, before the software is used to find solutions for the problem at hand. With so many research papers published on different ways of formulating plate and shell elements, it would be impossible to include all variations in this book; only the most fundamental ones related to plate element are listed here. Earlier in this chapter, a second order, 4-term polynomial equation is used to interpolate the in-plane displacements u and v for a plane stress element. The selection of the interpolation equation (a1 þ a2x þ a3y þ a4xy) makes this element type a bilinear element, as explained previously. In addition to the xy term, there are other second order terms: x2 and y2 shown in Fig. 2.22 (right). The triangular shaped figure, shown in Fig. 2.22 (left), is commonly known as Pascal’s triangle, which is credited to Blaise Pascal (Jun. 1623eAug. 1662), who published the concept. It is worth noting that in a number of old civilizations, such as China, India, and Persia, the same concept was reported long before Pascal’s publication. In addition to the triangle, Blaise Pascal was noted for his discoveries that the sum of the three internal angles of a triangle is equal to 180 degrees and the atmospheric pressure decreases as the height increases, because a vacuum exists in outer space. Pascal’s triangle is used to determine the number of terms in a complete nth order, two variables polynomial and the coefficients of all components with the form of (x þ y)n. First, the Pascal’s triangle on the right-hand side of Fig. 2.22 can be used to determine the number of complete, two-variable polynomial terms when given the order of the equation. For instance, a complete polynomial with the exponent n ¼ 0 has only one constant term (a0), which corresponds to the first row in the triangle. For n ¼ 1, there are three terms (a0 þ a1x þ a2y), corresponding to the first and second rows in the triangle. For a second order exponent (n ¼ 2), the complete Exponent 0 1

Pascal's Triangle 1 1 1

2

1

3

1

4 5

1 1

2 3

4 5

X 1

3 6

10

(X+Y) 1 X 1 4

10

X 1

5

X 1

X

Y XY

XY XY

XY

FIGURE 2.22 Pascal’s triangle and corresponding polynomial terms for (x þ y)n.

Y XY

XY XY

Y XY

XY

Y XY

Y

97

98

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

polynomial has the form of a0 þ a1x þ a2y þ a3x2þa4xy þ a5y2, as seen in the first three rows. This same procedure can be used to determine the complete polynomial terms for any two-variable equations with higher order exponents. This feature is used to select the appropriate terms to form an interpolation function needed to derive the element shape functions. For instance, an 8-node serendipity plane element requires eight polynomial terms and the eight associated coefficients needed to calculate the eight shape functions. However, the complete polynomial for a second order exponent (n ¼ 2) has a total of six terms, which is two terms short for a serendipity element. On the other hand, a third order exponent (n ¼ 3) has a total of ten terms, which is two terms more than that needed. If all six terms in the second order exponent are selected, code developers can decide which two terms selected from (x3, x2y, xy2, y3) are used to form the interpolation function. Alternatively, the code developer may want to choose combinational forms (e.g., x3 þ x2y and xy2 þ y3) for use in the calculation of the element shape functions. As such, shape functions for a serendipity element may differ a great deal due to different choices of polynomial terms. The second application of Pascal’s triangle is to determine the coefficients (a0, a1, a2, ., an) for a partial polynomial. A partial polynomial is related to only the terms within the same order. For example, a second order two-variable polynomial (x þ y)2 is related to only the x2, xy and y2 terms (i.e., the third row of the right Pascal’s triangle). The coefficients for the exponent n ¼ 2 are found in the third row of the left Pascal’s triangle, where the values are 1, 2, and 1. Hence, the resulting second order partial polynomial is (x þ y)2 ¼ x2þ2xy þ y2. Eq. (2.96) illustrates this concept up to the third order. If the component polynomials include coefficients, such as (2x þ y)2, the corresponding coefficients can be easily obtained by replacing x with 2x in Eq. (2.96). Based on this concept, we find (2x þ y)2 ¼ (2x)2 þ 2(2x) y þ y2 ¼ 4x2 þ 4xy þ y2. ðx þ yÞ0 ¼ 1; ðx þ yÞ1 ¼ x þ y; ðx þ yÞ2 ¼ x2 þ 2xy þ y2 ; ðx þ yÞ3 ¼ x3 þ 3x2 y þ 3xy2 þ y3 ; etc.

(2.96)

2.5.4.2 Select Polynomial Functions to Interpolate a Four-Node Plate Element For a 4-node plate element, there are a total of 12 DOFs, because each node has three DOFs. As we can see from the previous example using Pascal’s triangle, there are only 10 xy polynomial terms for a complete third order exponent. To obtain the 12 polynomial terms needed for the 4-node plate element, two additional terms from the fourth order polynomial must be involved. Because the complete fourth order polynomial has 15 xy terms, educated selections must be made based on individual or combined xy polynomial terms to obtain the necessary 12 terms. As mentioned earlier, different choices of polynomial terms will yield different interpolation equations and hence different shape functions. In this section, only the selection listed below is discussed. If you are solving problems in which plate and shell

2.5 Element Shape Functions and [B] Matrix

elements are of great importance, you should consider reading more research papers regarding the pros and cons of different polynomial terms selected to make up the interpolation function. For now, assume that a 12-term polynomial used to interpolate the vertical deflection w(x,y) of a 4-node plate element has the form wðx; yÞ ¼ a1 þ a2 x þ a3 y þ a4 x2 þ a5 xy þ a6 y2 þ a7 x3 þ a8 x2 y þ a9 xy2 þ a10 y3 þ a11 x3 y þ a12 xy3 . (2.97) Then, equations for the two slopes (in-plane rotations) are vwðx; yÞ ¼ a3 þ a5 x þ 2a6 y þ a8 x2 þ 2a9 xy þ 3a10 y2 þ a11 x3 þ 3a12 xy2 ; and vy (2.98) vwðx; yÞ ¼ a2 þ 2a4 x þ a5 y þ 3a7 x2 þ 2a8 xy þ a9 y2 þ 3a11 x2 y þ a12 y3 . (2.99) vx Written in matrix form, these three equations are fwðx; yÞg31 ¼ ½F312 fag121 0

8 9 > > wðx; yÞ > > > > > > 2 > > > > > > 1 x y x2 xy > > > < vwðx; yÞ > = 6 6 vy > ¼ 6 0 0 1 0 x > 4 > > > > > > > > > > 0 1 0 2x y > vwðx; yÞ > > > > : vx > ;

y2

x3

x2 y

xy2

y3

2y

0

x2

2xy

3y2

0

3x2

2xy

y2

0

9 a1 > > > > > > > > a2 > > > > > > > > a3 > > > > > > > > > > a4 > > > > > > a5 > > > > > > > > > a6 > = ; > a7 > > > > > > > > > a8 > > > > > > > > a9 > > > > > > > > > a10 > > > > > > > > a11 > > > > > > > ; a12 (2.100)

8 > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > 3 3 3 > > x y xy > > < 7> 7 x3 3xy2 7 5> > > > 2 3 > > 3x y y > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > :

where {w(x,y)}31 represents the generalized displacements (a vertical deflection and two in-plane rotations) at any point (x,y) within the element.

99

100

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

2.5.4.3 Identify 12 Constants for the Interpolation Polynomial We begin to find the 12 constants a1 through a12 by selecting the origin of the coordinate system to coincide with P1 and letting P2 be located at a distance on the positive x-axis, as shown in Fig. 2.23. To determine these 12 constants, let us first select the edge P1eP2. As we can see, all y-coordinates on this edge are zero. Hence, Eq. (2.97) becomes wðx; yÞ ¼ a1 þ a2 x þ a4 x2 þ a7 x3 .

(2.101)

We now consider a 4-node plate element as an analogue to a 2-node beam element. For the 2-node, 1D beam element described in Section 2.5.1, each node has two DOFs: vertical deflection w and rotation qy. For a 4-node, 2D plate element located on the xey plane, each node allows three DOFs: vertical deflection along the z-axis and two in-plane rotations, qy about the y-axis and qx about the x-axis. Eq. (2.101) describes the vertical deflection on the P1eP2 edge of a 4-node plate; hence, it is analogous to a beam element. Next, we take the derivative of the vertical deflection w, shown in Eq. (2.101), with respect to x at the same edge P1eP2 (where y ¼ 0). vwðx; yÞ ¼ qy ¼ a2 þ 2a4 x þ 3a7 x2 (2.102) vx Eqs. (2.101) and (2.102) describe w and qy of a 4-node plate element on the edge of y ¼ 0. These two equations are identical to the Hermite interpolation, shown in Eq. (2.25), used to formulate a 2-node beam element. Thus, the four needed constants, a1, a2, a4, and a7, can be solved as previously described, by applying the known nodal deflections (w1 and w2) and the slopes for rotation about the y-axis (q1y and q2y) to the interpolation function. As a quick summary, we have used two vertical deflections (w1 and w2) and two rotations about the y-axis (q1y and q2y) to find the four constants a1, a2, a4,

FIGURE 2.23 A plate of thickness h is represented by a 4-node surface element situated on the xey plane. Each node allows three DOFs: the vertical deflection wi, rotation qix about the xaxis, and rotation qiy about the y-axis.

2.5 Element Shape Functions and [B] Matrix

and a7. Let us now consider the slope vw vy , as shown in Eq. (2.98), at the same edge (P1eP2), where y ¼ 0. This results in vwðx; yÞ ¼ qx ¼ a3 þ a5 x þ a8 x2 þ a11 x3 . vy

(2.103)

We have four unknowns in Eq. (2.103), but only two slopes for the rotation about the x-axis (q1x and q2x) are left to be used. Because there are more unknowns than available equations, the slopes about the x-axis rotation have a potential for discontinuity. Therefore, we can conclude that the selection of the polynomial terms described in Eq. (2.97) allows a deflection continuity along the interelement boundaries, but not necessarily a continuous slope across the adjacent elements. As such, this formulation method may not provide the most accurate solutions. Despite this deficiency, it has been shown that plate elements formed from this set of polynomial terms provide acceptable solutions in a great number of cases. From Fig. 2.23 we can see that a positive vertical deflection (w) will generate a rotation about the “negative” y-axis. In other words, a positive w results in a counterclockwise rotation about the y-axis, which is opposite to the clockwise rotation shown based on the right-hand rule to represent “positive” rotation. We do not need to add any sign to the rotation about the x-axis because a positive w generates a rotation in the same way as that designated by the right-hand rule. As such, a negavw ðx;yÞ

tive sign is usually added to the matrix of the vx portion (i.e., the third row) of Eq. (2.100) to have matching sign conventions. Note that Eq. (2.100) is a 312 matrix that can be used to find w, qx, and qy at any point (x, y) within the plate element. Because we have four nodal points (x1, y1), (x2, y2), (x3, y3), and (x4, y4), inserting these four sets of (xi, yi) values into Eq. (2.100) yields a 1212 matrix. We then have 12 equations (1212 matrix) to solve for 12 unknowns, a1 through a12. We cannot identify a special set of nodal points that allows us to express the 12 constants in an elegant way, such as that shown for the shape functions of a 4-node rectangular element. Hence, sections below describe a numerical method to calculate the 12 constants. An alternative method for finding the 12 constants is to first solve the Hermite interpolation to identify a1, a2, a4, and a7, and then reduce the matrix size to 88 before finding the remaining eight constants from eight equations. Once all constants are found, they can be inserted back into Eq. (2.97), where the terms are arranged by nodal values {w1, q1x, q1y, w2, q2x,., q4y}, before the element shape functions are determined. Even with the number of equations being reduced to eight for finding the eight unknowns, manual calculations would be quite lengthy for obtaining this set of shape functions for a plate element. Use of a computer program is advised for finding the solutions. Procedures needed to write this computer program are outlined below. After finding the 12 constants, the element shape functions can be obtained.

101

102

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

2.5.4.4 Find Shape Functions for a 4-Node Plate Element Step 1: vw From Eq. (2.100), we can identify w; vw vy ; and vx (a vertical deflection and two rotations) as functions of 12 constants for each of the 4 nodes, because the nodal coordinates are known physical values. At node 1, the nodal coordinate is assumed to be (x1, y1), and the corresponding generalized displacements are (w1, q1x, q1y). Copying from Eq. (2.100), the top three rows of Eq. (2.104) can be written. As described earlier, we added a negative sign to the third row of this matrix to have a consistent way to enable the same direction of rotation due to the same sign of deflection. By repeating the same procedures for node 2 (x2, y2), we have the next three equations for (w2, q2x, q2y). Similarly, we obtain the six additional equations for (w3, q3x, q3y) and (w4, q4x, q4y). Eq. (2.104) shows the combination of these 12 equations in matrix form. In other words, the 312 matrix shown in Eq. (2.100) will become a 1212 matrix, and this new matrix is assigned as the [G] matrix.

nodalfwg121

2

1

6 6 60 6 6 60 6 6 6 61 6 ¼6 60 6 6 6 6: 6 6 6 6: 4 0

8 9 w1 > > > > > > > > > > > > > > q > 1x > > > > > > > > > > > > > > > q 1y > > > > > > > > > > > < w2 > = ¼ ½G1212 fag121 0 > > > > > > q2x > > > > > > > > > > > > : > > > > > > > > > > > > > > : > > > > > > > > > > : ; q4y

x1

y1

x1 2

x1 y1

y1 2

x1 3

x1 2 y1

x1 y1 2

y1 3

x1 3 y1

0

1

0

x1

2y1

0

x1 2

2x1 y1

3y1 2

x1 3

1

0

2x1

y1

0

3x1 2

2x1 y1

y1 2

0

3x1 2 y1

x2

y2

x2 2

x2 y2

y2 2

x2 3

x2 2 y2

x2 y2 2

y2 3

x2 3 y2

0

1

0

x2

2y2

0

x2 2

2x2 y2

3y2 2

x2 3

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

1

0

2x4

y4

0

3x4 2

2x4 y4

y4 2

0

3x4 2 y4

38 9 a1 > > > > > > 7> > > > > > 27 > a2 > > > 3x1 y1 7 > > 7> > > 7> > > > > 3 7> > > a y1 7> 3 > > > > 7> > > > > > > 3 7 = < a4 > x2 y2 7 7 7 >a > > 3x2 y2 2 7 > 7> > > 5 > > 7> > > > 7> > > > > 7 : : 7> > > > > > > 7> > > > 7> > : > > : 7> > > > > 5> > > > ; : 3 a12 y4 (2.104) x1 y1 3

2.5 Element Shape Functions and [B] Matrix

Step 2: We can easily find [G]1 numerically from Eq. (2.104) because (x1, y1) through (x4, y4) are known. By multiplying [G]1 to both sides of Eq. (2.104), we find the 12 constants, a1 through a12, as shown in Eq. (2.105). fag ¼ ½G1 f w1

q1x

q1y

.

q4y gT

(2.105)

Step 3: We now write the generalized displacements w at any point (x, y) within the element as a function of nodal displacements by inserting Eq. (2.105) into Eq. (2.100), as shown in Eq. (2.106). 8 9 > > wðx; yÞ > > > > > > > > > > > < vwðx; yÞ > = ¼ ½F312 fag121 fwðx; yÞg ¼ > vy > > > > > > > vwðx; yÞ > > > > > : vx > ; ¼ ½F312 ½G1 1212 f w1

q1x

q1y

.

q4y gT 121

(2.106)

By definition, the element shape functions are used to derive a physical value at any point within the element; that is, fwðx; yÞg ¼ ½Nf w1 q1x q1y . q4y gT . Thus, the shape functions for a 4-node plate element, with the interpolation function selected as described above, is expressed as [N] ¼ [F][G]1.

2.5.4.5 Determine StraineDisplacement Matrix As described for a beam element, the strainedisplacement relationship for a plate element is equivalent to the curvatureedisplacement relationship. By recalling how the curvature was derived for a beam element, we can deduce that the generalized curvatures for a plate element have the form v2 w ¼ 2a4 þ 6a7 x þ 2a8 y þ 6a11 xy; vx2

(2.107)

v2 w ¼ 2a6 þ 2a9 x þ 6a10 y þ 6a12 xy; and vy2

(2.108)

2v2 w ¼ 2a5 þ 4a8x þ 4a9 y þ 6a11 x2 þ 6a12 y2 . vxvy

(2.109)

kx ¼ ky ¼ kxy ¼

Note that in some textbooks, a negative sign is assigned to the curvatures, but for consistency with the definition used in the formulation of the [B] matrix for the beam element, the negative sign convention is not adopted here.

103

104

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

Step 4: Rewrite the generalized curvatures shown in Eqs. (2.107)e(2.109) in matrix form as a function of the 12 constants, a1 through a12. This new matrix, [H], consists of all the coefficients, as shown in Eq. (2.110). 8 9 > < kx > = fkg ¼ ½Hfag0 ky > > : ; kxy

2

0 6 6 ¼6 0 4 0

2

0

0

2

0

0

6x 2y

0

0

6xy

0

0

0

0

2

0

0

2x

6y

0

0

0

0

2

0

0

4x

4y

0

6x2

3 a1 6a 7 6 2 7 7 36 6 a3 7 0 6 7 76 a 7 76 4 7 6xy 76 7 56 : 7 6 7 7 6y2 6 6 : 7 6 7 4 : 5 a12 (2.110)

Step 5: Insert {a} from Eq. (2.105) into Eq. (2.110), to obtain Eq. (2.111). fkg ¼ ½Hfag ¼ ½H½G1 f w1 ¼ ½Bf w1

q1x

q1y

q1x

q1y

. q4y gT

. q4y gT (2.111)

Hence, the strainedisplacement or curvatureedisplacement matrix [B] for a 4node plate element can be written as [H][G]1.

2.5.5 3D, 4-NODE SHELL ELEMENT A 4-node shell element is usually formed by the superposition of a 4-node plane stress element and a 4-node plate element. Because a plane stress element has two DOFs per node and a plate element has three DOFs per node, a shell element has a total of five DOFs per node. The only DOF not accounted for in a shell element is the rotation about the z-axis (i.e., in-plane rotation), also known as the drilling DOF. This five-DOFs shell element is problematic, because a real-world structure usually involves connecting a shell to a beam or another shell not in the same reference surface. Hence, it is more convenient to artificially make a shell element having a total of six DOFs (three translational and three rotational DOFs). This treatment equips a 4-node shell element with a total of 24 DOFs per element. So, the order of the stiffness matrix for this augmented shell element is 24  24. Eq. (2.112)

2.5 Element Shape Functions and [B] Matrix

shows the results of superposition of a plane stress element, a plate element, and zero stiffness for the four drilling DOFs. 3 2 0 0 ½kplane 88 7 6 (2.112) ½k2424 ¼ 4 0 ½kplate 1212 0 5 0

0

044

The shell element formulated as shown above contains no stiffness for in-plane rotation, hence Eq. (2.112) is a singular matrix. To resolve this problem, a small amount of artificial torsional stiffness is supplemented in some software packages. The obvious question is how much artificial stiffness is small enough to ensure accuracy from the model. Because there may be different polynomial terms selected for a plate element during formulation, there are several plate elements that exist. Similarly, there are numerous shell elements available. One such method is presented by Rengarajan et al. (1995). However, discussion about different ways a shell element can be formulated is beyond the scope of this book. Advanced readers may want to look into scientific literature or study the theoretical manual provided by the software vendor in order to better understand the pros and cons regarding how to select the proper element type.

2.5.6 3D, 8-NODE TRILINEAR ELEMENT SHAPE FUNCTIONS Fig. 2.24 shows an 8-node solid brick element with dimensions of 2a by 2b by 2c. The origin of this brick element is located at its geometric center. Again, the numbering sequence for this element is important. In this case, P1 is located at the bottom layer, at coordinates (a, b, c), while P2, P3, and P4 are formed through a counterclockwise rotation. For the top layer, P5 is located at coordinates (a, b, c), and P6 through P8 are formed using the same counterclockwise rotation. Assume the interpolation function for this solid element has the form 4ðx; y; zÞ ¼ a1 þ a2 x þ a3 y þ a4 z þ a5 xy þ a6 yz þ a7 zx þ a8 xyz ¼

8 X

Ni 4i .

i¼1

(2.113) This is a trilinear interpolation equation, because it varies linearly along the x-axis if both the y and z values are held constant. Similarly, it behaves linearly along the y-axis if the x and z values are held constant, and linearly along the z-axis if the x and y values are held constant. By now, we are familiar with a couple of methods used to identify the element shape functions. Although it is quite tedious, we have enough information thus far to be able to derive this set of shape functions. As such, detailed

105

106

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

FIGURE 2.24 An 8-node trilinear element of size 2a  2b  2c. The origin of the local coordinate system is located at the geometric center. The nodes are arranged such that P1 through P4 are at the bottom layer while P5 through P8 are at the top layer. In both layers, these nodes are arranged in a counterclockwise manner.

derivations are not provided. Instead, only the resulting element shape functions are listed below. ða  xÞðb  yÞðc  zÞ ða þ xÞðb  yÞðc  zÞ ; N2 ¼ 8abc 8abc ða þ xÞðb þ yÞðc  zÞ ða  xÞðb þ yÞðc  zÞ N3 ¼ ; N4 ¼ 8abc 8abc (2.114) ða  xÞðb  yÞðc þ zÞ ða þ xÞðb  yÞðc þ zÞ N5 ¼ ; N6 ¼ 8abc 8abc ða þ xÞðb þ yÞðc þ zÞ ða  xÞðb þ yÞðc þ zÞ N7 ¼ ; N8 ¼ 8abc 8abc We can use the same shape functions for all three translational DOFs and rewrite Eq. (2.113) in matrix format as 8 9 u1 > > > > > > > > > > v1 > > > > 2 3 8 9 > > > N1 0 0 N2 0 0 : : 0 > w1 > > uðx; y; zÞ > > > > < = < = 6 7 6 7 vðx; y; zÞ ¼ 4 0 N1 0 0 N2 0 : : 0 5 u2 . > > > > > > : ; > > wðx; y; zÞ > : > > 0 0 N1 0 0 N2 : : N8 324> > > > > : > > > > > > > : > ; w8 241 (2.115) N1 ¼

2.5 Element Shape Functions and [B] Matrix

107

By inserting this equation into Eq. (1.11), the strainedisplacement equation, the following 3D strainedisplacement equations can be established: 8 9 8 9 εxx > > εxx > > > > > > > > > > > > > > > > > > > > > > ε ε yy yy > > > > > > > > > > > > > > > < ε = < ε > = zz zz ¼ > > gxy > 2εxy > > > > > > > > > > > > > > > > > > > > > 2εyz > > gyz > > > > > > > > > > > > > : : ; > ; gzx 2εzx 2 3 v 0 0 6 vx 7 6 7 8 9 6 7 u1 > 6 0 v 0 7 > > vy 6 7 > > > > 6 7 > > v1 > > > 6 7 2 > > 3 > > > > 6 7 v N 0 0 N 0 0 : : 0 > > 1 2 w > 6 0 0 vz 7 1> < = 6 7 6 7 6 7 6 7 ¼ 6 v v 0 7 4 0 N1 0 0 N2 0 : : 0 5 u2 . > > > > 6 vy vx 7 > > : > 6 7 0 0 N1 0 0 N2 : : N8 324> > > > > 6 7 > > > > 6 0 v v 7 > > : > > 6 7 vz vy > : > ; 6 7 w8 241 6 7 6 7 4 v v 5 0 vx vz 63 (2.116) Finally, the strainedisplacement matrix [B] can be written as 3 2 v 0 0 7 6 vx 7 6 7 6 6 0 v 0 7 vy 7 6 7 6 7 2 6 7 6 v N1 0 0 N2 0 0 : : 6 0 0 vz 7 7 6 6 7 6 0 N1 0 ½B ¼ 6 0 N2 0 : : v v 7 4 6 vy vx 0 7 6 7 6 0 0 N1 0 0 N2 : : 7 6 6 0 v v 7 6 vz vy 7 7 6 7 6 7 6 5 4 v v 0 vx vz 63

0

3

7 0 7 . 5 N8 324

(2.117)

108

CHAPTER 2 Meshing, Element Types, and Element Shape Functions

EXERCISES 1. Consider a 4-node surface element with P1 (0, 0), P2 (2, 0), P3 (2, 3.1), and P4 (0, 3). Find the element shape functions in this global coordinate system. 2. Consider a 4-node bilinear element having nodal coordinates P1 (0, 0), P2 (6, 0), P3 (6, 4), and P4 (0, 4). (i) Calculate the average shear strain after deformation which causes the coordinates of P3 and P4 to become (6.1, 4.0) and (0.1, 4.0), (ii) Assuming the origin of the coordinate system is located at P1, show step-by-step derivations of the four element shape functions [N], (iii) Determine matrix [B] based on the equation 2 the strainedisplacement 2 3 3 8 9 u1 > v v > > > 9 6 6 > > 07 0 78 P 4 > 8 9 6 vx > 6 vx > 7 7> v1 > > > > > > > > 6 6 7 7 N u i i > > > > >   6 < < < εxx > = = 6 7 7 v7 u v7 1 u2 = 6 6 0 εyy ¼ 6 0 , and ¼ ½B ¼6 7 7 X 4 vy 7 v vy 7> > > > > 6 : > > > : > ; 6 > > > > 6 6 7 7> gxy N v > : > ; i i > 6v v7 6v v7 > > : > > > 4 4 5 5 1 > : > ; vy vx vy vx v4 (iv). Calculate the element strain from the [B] matrix when the coordinates of P3 and P4 become (6.1, 4.0) and (0.1, 4.0). 3. Create a table of all the 1D and 2D elements, and some 3D elements with the number of degrees of freedom, type of degrees of freedom, what kinds of loads are resisted, and the stiffness matrix dimensions. 4. Find the shape functions for a 1D bar element with the two points P1 and P2 located at coordinates x ¼ 3 and x ¼ 5. 5. Find the strainedisplacement [B] matrix of the element in problem 4. 6. Find the shape functions of a beam element that has endpoints of x ¼ 3.5 and x ¼ 7.2. 7. Create graphs of the four shape functions of a beam element that has endpoints of x ¼ 3.5 and x ¼ 7.2. 8. Find the [B] matrix of the beam element used in problems 6 and 7. 9. If a constant strain triangle element has points P1 (0, 0), P2 (4, 0), and P3 (2, 5), what are the shape functions of this element? 10. Find a polynomial equation (x) that passes through three points on the xey plane with coordinates of (0.6, 12), (6.2, 3.9), and (10.1, 9).

References

11. A room 15  40 m in size has points P1 through P4 arranged in a counterclockwise manner in the corners of the room. The air pressures at these points are 100, 105, 108, and 98, respectively. Find and plot the contour lines of 100, 103, and 106 kPa over the room area. 12. Use Lagrange interpolation to find the shape functions for a 3D, 8-node trilinear element with dimensions 2a  2b  2c.

REFERENCES Hallquist, J.O., 2006. LS-DYNA Theoretical Manual. ISBN:0-9778540-0-0, Available through LSTC website: http://lstc.com/download/manuals. Jeffreys, H., Jeffreys, B.S., 1988. Lagrange’s interpolation formula. x9.011. In: Methods of Mathematical Physics, third ed. Cambridge University Press, Cambridge, England, p. 260. Love, A.E.H., 1888. The small free vibrations and deformation of a thin elastic shell. Philosophical Transactions of the Royal Society of London A 179, 491e546. Mindlin, R.D., 1951. Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates. ASME Journal of Applied Mechanics 18, 31e38. O’Connor, J.J., Robertson, E.F., 2017. MacTutor History of Mathematics Archive. URL: http://www-history.mcs.st-andrews.ac.uk/index.html. Reissner, E., 1945. The effect of transverse shear deformation on the bending of elastic plates. ASME Journal of Applied Mechanics 12, A68eA77. Rengarajan, G., Aminpour, M.A., Knight Jr., N.F., 1995. Improved assumed-stress hybrid shell element with drilling degrees of freedom for linear stress, buckling and free vibration analyses. International Journal for Numerical Methods in Engineering 38, 1917e1943. Schneiders, R., 2017. URL: http://www.robertschneiders.de/meshgeneration/software.html. Wolfram Math World, 2017. URL: http://mathworld.wolfram.com.

109