Introduction
In this topic, we develop, in an easy to understand manner, finite element equations for the stress analysis of two-dimensional (2D) solids subjected to external loads.The element developed is called a 2D solid element that is used for structural problems where the loading-and hence the deformation-occur within a plane. Though no real life structure can be truly 2D, experienced analysts can often idealize many practical problems to 2D problems to obtain satisfactory results by carrying out analyses using 2D models, which can be very much more efficient and cost-effective compared to conducting full 3D analyses. In engineering applications, there are ample practical problems that can be modelled as 2D problems.For example, if we have a plate structure with loading acting in the plane of the plate as in Figure 7.1, we need to use 2D plane stress elements. When we want to model the effects of water pressure on a dam, as shown in Figure 7.2, we have to use 2D plane strain elements.
Figure 7.1. Atypical 2D plane stress problem.
Figure 7.2. Atypical 2D plane strain problem.
Note that in Figure 7.1, plane stress conditions are usually applied to structures that have a relatively small thickness as compared to its other dimensions. Due to the absence of any off-plane external force, the normal stresses are negligible, which leads to a plane stress situation. In cases where plane strain conditions are applied, as in Figure 7.2, the thickness of the structure (in the z direction) is relatively large as compared to its other dimensions, and the loading (pressure) is uniform along the elongated direction. The deformation is, therefore, approximated to be the same throughout its thickness.The formulations for plane stress and plane strain problems are very much the same, except for the difference in the material constant matrix.
A 2D solid element, be it plane strain or plane stress, can be triangular, rectangular or quadrilateral in shape with straight or curved edges. The most often used elements in engineering practice are linear. Quadratic elements are also used for situations that require high accuracy in stress, but they are less often used for practical problems. Higher order elements have also been developed, but they are less often used except for certain specific problems. The order of the 2D element is determined by the order of the shape functions used. A linear element uses linear shape functions, and therefore the edges of the element are straight. A quadratic element uses quadratic shape functions, and their edges can be curved. The same can be said for elements of third order or higher.
In a 2D model, the elements can only deform in the plane where the model is defined, and in most situations, this is the x-y plane. At any point, the variable, that is the displacement, has two components in the x and y directions, and so do the external forces. For plane strain problems, the thickness of the true structure is usually not important, and is normally treated as a unit quantity uniformly throughout the 2D model. However, for plane stress problems, the thickness is an important parameter for computing the stiffness matrix and stresses. Throughout this topic, it is assumed that the elements have a uniform thickness of A. If the structure to be modelled has a varying thickness, the structure needs to be divided into small elements, where in each element a uniform thickness can be used. On the other hand, formulation of 2D elements with varying thicknesses can also be done easily, as the procedure is similar to that of a uniform element. Very few commercially available software packages provide elements of varying thickness.
The equation system for a 2D element will be more complex as compared with the 1D element because of the higher dimension. The procedure for developing these equations is, however, very similar to that for the 1D truss elements.These steps can be summarized in the following three-step procedure:
1. Construction of shape functions matrix N that satisfies Eqs. (3.34) and (3.41).
2. Formulation of the strain matrix B.
3. Calculation of ke, me, and fe using N and B and Eqs. (3.71), (3.75) and (3.81).
We shall be focusing on the formulation of three types of simple but very important elements: linear triangular, bilinear rectangular, and isoparametric linear quadrilateral elements. Once the formulation of these three types of element is understood, the development of other types of elements of higher orders is straightforward, because the same techniques can be utilized. Development of higher order elements will be discussed at the end of this topic.
Linear Triangular Elements
The linear triangular element was the first type of element developed for 2D solids. The formulation is also the simplest among all the 2D solid elements. It has been found that the linear triangular element is less accurate compared to linear quadrilateral elements. For this reason, it is often thought to be ideal to use quadrilateral elements, but the reality is that the triangular element is still a very useful element for its adaptivity to complex geometry. Triangular elements are normally used when we want to mesh a 2D model involving complex geometry with acute corners. In addition, the triangular configuration with the simplest topological feature makes it easier to develop meshing processors. Nowadays, analysts are hoping to use a fully automated mesh generator to perform the complex task of analysis that needs repeated or even adaptive re-meshing. Most automated mesh generators can only create triangular elements. There are automated mesh generators that can generate a quadrilateral mesh, but they still use triangular elements as patches for difficult situations, and end up with a mesh of mixed elements. Hence, whether we like it or not, we still have to use triangular elements for many practical engineering problems.
Consider a 2D model in the x-y plane, shown schematically in Figure 7.3. The 2D domain is divided in a proper manner into a number of triangular elements.In a mesh of linear triangular elements, each triangular element has three nodes and three straight edges.
Figure 7.3. Rectangular domain meshed with triangular elements.
Figure 7.4. Linear triangular element.
Field Variable Interpolation
Consider now a triangular element of thickness h. The nodes of the element are numbered 1, 2 and 3 counter-clockwise, as shown in Figure 7.4. For 2D solid elements, the field variable is the displacement, which has two components (u and v), and hence each node has two Degrees Of Freedom (DOFs). Since a linear triangular element has three nodes, the total number of DOFs of a linear triangular element is six. For the triangular element, the local coordinate of each element can be taken as the same as the global coordinate, since there are no advantages in specifying a different local coordinate system for each element.
Now, let us examine how a triangular element can be formulated. The displacement U is generally a function of the coordinates x and y, and we express the displacement at any point in the element using the displacements at the nodes and shape functions. It is therefore assumed that (see Section 3.4.2)
where the superscript h indicates that the displacement is approximated, and de is a vector of the nodal displacements arranged in the order of
and the matrix of shape functions N is arranged as
in which Ni (i = I, 2, 3) are three shape functions corresponding to the three nodes of the triangular element. Equation (7.1) can be explicitly expressed as
which implies that each of the displacement components at any point in the element is approximated by an interpolation from the nodal displacements using the shape functions. This is because the two displacement components are basically independent from each other. The question now is how can we construct these shape functions for our triangular element that satisfies the sufficient requirements: delta function property; partitions of unity; and linear field reproduction.
Shape Function Construction
Development of the shape functions is normally the first, and most important, step in developing finite element equations for any type of element. In determining the shape functions Ni (i = 1, 2, 3) for the triangular element, we can of course follow exactly the standard procedure described in Sections 3.4.3 and 4.2.1, by starting with an assumption of the displacements using polynomial basis functions with unknown constants. These unknown constants are then determined using the nodal displacements at the nodes of the element. This standard procedure works in principle for the development of any type of element, but may not be the most convenient method. We demonstrate here another slightly different approach for constructing shape functions. We start with an assumption of shape functions directly using polynomial basis functions with unknown constants. These unknown constants are then determined using the property of the shape functions. The only difference here is that we assume directly the shape function instead of the displacements. For a linear triangular element, we assume that the shape functions are linear functions of x and y. They should, therefore, have the form of
where ai, bi and Ci (i = 1, 2, 3) are constants to be determined. Equation (7.5) can be written in a concise form,
We write the shape functions in the following matrix form:
where a is the vector of the three unknown constants, and p is the vector of polynomial basis functions (or monomials). Using Eq. (3.21), the moment matrix P corresponding to basis p can be given by
Note that the above equation is written for the shape functions, and not for the displacements. For this particular problem, we use up to the first order of polynomial basis. Depending upon the problem, we could use higher order of polynomial basis functions. The complete order of polynomial basis functions in two-dimensional space up to the «th order can be given by using the so-called Pascal triangle, shown in Figure 3.2. The number of terms used in p depends upon the number of nodes the 2D element has. We usually try to use terms of lowest orders to make the basis as complete as possible in order. It is also possible to choose specific terms of higher orders for different types of elements. For our triangular element there are three nodes, and therefore the lowest terms with complete first order are used, as shown in Eq. (7.9). The assumption of Eq. (7.5) implies that the displacement is assumed to vary linearly in the element. In Eq. (7.8) there is a total of nine constants to be determined. Our task now is to determine these constants.
Therefore, we can expect that the complete linear basis functions used in Eq. (7.9) guarantee that the shape functions to be constructed satisfy the sufficient requirements for FEM shape functions. What we need to do now is simply impose the delta function property on the assumed shape functions to determine the unknown constants ai, bi and Ci.
The delta functions property states that the shape function must be a unit at its home node, and zero at its remote nodes. For a two-dimensional problem, it can be expressed as
For a triangular element, this condition can be expressed explicitly for all three shape functions in the following equations. For shape function Ni, we have
This is because node 1 at (x\, yO is the home node of N1, and nodes 2 at (x2, y2) and 3 at (x3, y3) are the remote nodes of N1. Using Eqs. (7.5) and (7.12), we have
Solving Eq. (7.13) for a1, b1 and C1, we obtain
where Ae is the area of the triangular element that can be calculated using the determinant of the moment matrix:
Note here that as long as the area of the triangular element is nonzero, or as long as the three nodes are not on the same line, the moment matrix P will be of full rank. Substituting Eq. (7.14) into Eq. (7.5), we obtain
which can be re-written as
This equation clearly shows that N1 is a plane in the space of (x, y, N) that passes through the line of 2-3, and vanishes at nodes 2 at (x2, y2) and3at(x3, y3). Thisplanealsopassesthe point of (x1,y1, 1) in space that guarantees the unity of the shape function at the home node. Since the shape function varies linearly within the element, N1 can then be easily plotted as in Figure 7.5(a). Making use of these features of N1, we can immediately write out the other two shape functions for nodes 2 and 3. For node 2, the conditions are
and the shape function N2 should pass through the line 3-1, which gives
which is plotted in Figure 7.5(b).
Figure 7.5. Linear triangular element and its shape functions.
For node 3, the conditions are
and the shape function N3 should pass through the line 1-2, and given by
which is plotted in Figure 7.5(c). Theprocess of determining these constants is basically simple, algebraic manipulation. The shape functions are summarized in the following concise form:
where the subscript i varies from 1 to 3, and j and k are determined by the cyclic permutation in the order of i, j, k. For example, if i = 1, then j = 2,k = 3. When i = 2, then j = 3, k = 1.
Area Coordinates
Another alternative and effective method for creating shape functions for triangular elements is to use so-called area coordinates L\,L2 and L3. The use of the area coordinates will immediately lead to the shape functions for triangular elements. However, we first need to define the area coordinates.
In defining Li, we consider a point Pat (x, y) inside the triangle, as shown in Figure 7.6, and form a sub-triangle of 2-3-P. The area of this sub-triangle is noted as Ai, and can be calculated using the formula
Figure 7.6. Definition of area coordinates.
The area coordinate L1 is then defined as
Similarly, for L2 we form sub-triangle 3-1-P with an area of A2 given by
The area coordinate L2 is then defined as
For L3, we naturally write
where A3 is the area of the sub-triangle 1-2-P, calculated using
It is very easy to confirm the unity property of the area coordinates L1,L2 and L3. First, they are partitions of unity, i.e.
that can be proven using the definition of the area coordinates:
Secondly, these area coordinates are of delta function properties. For example, L1 will definitely be zero if P is at the remote nodes 2 and 3, and it will be a unit if P is at its home node 1. The same arguments are also valid for L2 and L3.
These two properties are exactly those defined for shape functions. Therefore, we immediately have
The previous equation can also be easily confirmed by comparing Eqs. (7.16) with (7.25), (7.19) with (7.27), and (7.21) with (7.28). The area coordinates are very convenient for constructing higher order shape functions for triangular elements.
Once the shape function matrix has been developed, one can write the displacement at any point in the element in terms of nodal displacements in the form of Eq. (7.1). The next step is to develop the strain matrix so that we can write the strain, and hence the stress, at any point in the element in terms of the nodal displacements. This will further lead to the element matrices.
Strain Matrix
Let us now move to the second step, which is to derive the strain matrix required for computing the stiffness matrix of the element.There are only three major stress components,
in a 2D solid, and the corresponding strains,
can be expressed as
or in a concise matrix form,
where L is called a differential operation matrix, and can be obtained simply by inspection
of Eq. (7.33):
Substituting Eq. (7.1) into Eq. (7.34), we have
where B is termed the strain matrix, which can be obtained by the following equation once the shape function is known:
Equation (7.36) implies that the strain is now expressed by the nodal displacement of the element using the strain matrix. Equations (7.36) and (7.37) are applicable for all types of 2D elements.
Using Eqs. (7.3), (7.22), (7.23) and Eq. (7.37), the strain matrix B for the linear triangular element can be easily obtained, to have the following simple form:
It can be clearly seen that the strain matrix B for a linear triangular element is a constant matrix. This implies that the strain within a linear triangular element is constant, and thus so is the stress. Therefore, the linear triangular elements are also referred to as constant strain elements or constant stress elements. In reality, stress or strain varies across the structure, hence using linear triangular elements with a coarse mesh will result in a rather inaccurate stress or strain distribution. We would need to have a fine mesh of linear triangular elements in order to show an appropriate variation of stress or strain across the structure.
Element Matrices
Having obtained the shape function and the strain matrix, the displacement and strain (hence the stress) can all be expressed in terms of the nodal displacements of the element. The element matrices, like the stiffness matrix ke, mass matrix me.
The element stiffness matrix ke for 2D solid elements can be obtained using Eq. (3.71):
Note that the material constant matrix c has been given by Eqs. (2.31) and (2.32), respectively, for the plane stress and plane strain problems. Since the strain matrix B is a constant matrix, as shown in Eq. (7.38), and the thickness of the element is assumed to be uniform, the integration in Eq. (7.39) can be carried out very easily, which leads to
The element mass matrix me can also be easily obtained by substituting the shape function matrix into Eq. (3.75):
For elements with uniform thickness and density, we can rewrite Eq. (7.41) as
The integration of all the terms in the mass matrix can be carried out by simply using a mathematical formula developed by Eisenberg and Malvern (1973):
where
is the area coordinates for triangular elements that is the same as the shape function, as we have seen in Section 7.2.2. The element mass matrix me is found to be
The nodal force vector for 2D solid elements can be obtained using Eqs. (3.78), (3.79) and (3.81). Suppose the element is loaded by a distributed force f on the edge 2-3 of the element, as shown in Figure 7.4; the nodal force vector becomes
If the load is uniformly distributed,are constants within the element, so the above equation becomes
where l2-3 is the length of the edge 2-3 of the element.
Once the element stiffness matrix ke, mass matrix me and nodal force vector fe have been obtained, the global finite element equation can be obtained by assembling the element matrices by summing up the contribution from all the adjacent elements at the shared nodes.
