FEM for Two-Dimensional Solids (Finite Element Method) Part 2

Linear Rectangular Elements

Triangular elements are usually not preferred by many analysts nowadays, unless there are difficulties with the meshing and re-meshing of models of complex geometry. The main reason is that the triangular elements are usually less accurate than rectangular or quadrilateral elements. As shown in the previous section, the strain matrix of the linear triangular elements is constant, accounting for the inaccuracy. With advances in meshing algorithms, many models of complex geometry with sharp corners or curved edges can be modelled using quadrilateral elements. For the rectangular element, the strain matrix is not a constant, as will be shown in this section. This will provide a more realistic presentation in strain, and hence the stress distribution, across the structure. The formulation of the equations for the rectangular elements is simpler compared to the triangular elements, because the shape functions can form very easily due to the regularity in the shape of the rectangular element. The simple three-step procedure is applicable, and will be shown in the following sections.

Shape Function Construction

Consider a 2D domain. The domain is discretized into a number of rectangular elements with four nodes and four straight edges, as shown in Figure 7.7. As always, we number the nodes in each element 1, 2, 3 and 4 in a counter-clockwise direction. Note also that, since each node has two DOFs, the total DOFs for a linear rectangular element would be eight.

The dimension of the element is defined here as 2a x 2b x h. A local natural coordinate system (ξ, η) with its origin located at the centre of the rectangular element is defined. The relationship between the physical coordinate (x, y) and the local natural coordinate system (ξ, η) is given by

Equation (7.47) defines a very simple coordinate mapping between physical and natural coordinate systems for rectangular elements as shown in Figure 7.8. Our formulation can now be based on the natural coordinate system. The use of natural coordinates will make the construction of the shape functions and evaluation of the matrix integrations very much easier. This kind of coordinate mapping technique is one of the most often used techniques in the FEM. It is extremely powerful when used for developing elements of complex shapes.

We perform the field variable interpolation and express the displacement within the element as an interpolation of the nodal displacements using shape functions. The displacement vector U is assumed to have the form

where the nodal displacement vector de is arranged in the form

and the matrix of shape functions has the form

 

Figure 7.7. Rectangular domain meshed by rectangular elements.

 

where the shape functions Ni (i = I, 2, 3,4) are the shape functions corresponding to the four nodes of the rectangular element.

In determining these shape functions Ni (i = I, 2, 3, 4), we can follow exactly the same steps used in Sections 4.2.1 or 7.2.2, by starting with an assumption of the displacement or shape functions using polynomial basis functions with unknown constants. These unknown constants are then determined using the displacements at the nodes of the element or the property of the shape functions. The only difference here is that we need to use four terms of monomials of basis functions. As we have seen in Section 7.2.2, the process is quite troublesome and lengthy. For many cases, one often constructs shape functions simply by some ‘shortcut’ methods. One of these is by inspection, and utilizing the properties of shape functions.

Due to the regularity of the square element in the natural coordinates, the shape functions in Eq. (7.50) can be written out directly as follows, without going through the detailed process that we described in the previous section for triangular elements:

 

Figure 7.8. Rectangular element and the coordinate systems. (a) Rectangular element in physical system, (b) square element in natural coordinate system.

It is very easy to confirm that all the shape functions given in Eq. (7.51) satisfy the delta function property of Eq. (3.34). For example, for N3 we have

The same examination of N1, N2 and N4 will confirm the same property.

It is also very easy to confirm that all the shape functions given in Eq. (7.51) satisfy the partition of unity property of Eq. (3.41):

 

Equation (7.53) should be called a bilinear shape function to be exact, as it varies linearly in both the ξ and η directions. It varies quadratically in any direction other than these two ξ and η directions. Denoting the natural coordinates of node j by,    the bilinear shape function Nj can be re-written in a concise form:

Strain Matrix

Using the same procedure as for the case of the triangular element, the strain matrix B would have the same form as in Eq. (7.37), that is

It is now clear that the strain matrix for a bilinear rectangular element is no longer a constant matrix. This implies that the strain, and hence the stress, within a linear rectangular element is not constant.

Element Matrices

Having obtained the shape function and the strain matrix B, the element stiffness matrix ke. Using first the relationship given in Eq. (7.47), we have

Substituting Eq. (7.56) into Eq. (7.39), we obtain

The material constant matrix c has been given by Eqs. (2.31) and (2.32), respectively, for plane stress and plane strain problems. Evaluation of the integral in Eq. (7.57) would not be as straightforward, since the strain matrix B is a function of ξ and η. It is still possible to obtain the closed form for the stiffness matrix by carrying out the integrals in Eq. (7.57) analytically. In practice, we often use a numerical integration scheme to evaluate the integral, and the commonly used Gauss integration scheme will be introduced here. The Gauss integration scheme is a very simple and efficient procedure that performs numerical integral, and it is briefly outlined here.

Gauss Integration

Consider first a one-dimensional integral. Using the Gauss integration scheme, the integral is evaluated simply by a summation of the integrand evaluated at m Gauss points multiplied by corresponding weight coefficients as follows:

The locations of the Gauss points and the weight coefficients have been found for different m, and are given in Table 7.1. In general, the use of more Gauss points will produce more accurate results for the integration. However, excessive use of Gauss points will increase the computational time and use up more computational resources, and it may not necessarily give better results. The appropriate number of Gauss points to be used depends upon the complexity of the integrand. It has been proven that the use of m Gauss points gives the exact results of a polynomial integrand of up to an order of n = 2m – 1. For example, if the integrand is a linear function (straight line), we have 2m – 1 = 1, which gives m = 1. This means that for a linear integrand, one Gauss point will be sufficient to give the exact result of the integration.

Table 7.1. Gauss integration points and weight coefficients

m			Accuracy n
1	0	2	1
2	— 1/V3, 1/V3	1,1	3
3	—VÔ6, 0,V06	5/9, 8/9, 5/9	5
4	—0.861136, —0.339981,	0.347855, 0.652145,	7
	0.339981, 0.861136	0.652145, 0.347855
5	—0.906180, —0.538469, 0,	0.236927, 0.478629, 0.568889,	9
	0.538469, 0.906180	0.478629, 0.236927
6	—0.932470, —0.661209, —0.238619,	0.171324, 0.360762, 0.467914,	11
	0.238619, 0.661209, 0.932470	0.467914, 0.360762, 0.171324

If the integrand is of a polynomial of a third order, we have 2m — I = 3, which gives m = 2. This means that for an integrand of a third order polynomial, the use of two Gauss points will be sufficient to give the exact result. The use of more than two points will still give the same results, but takes more computation time. For two-dimensional integrations, the Gauss integration is sampled in two directions, as follows:

Figure 7.9(b) shows the locations of four Gauss points used for integration in a square region.

The element stiffness matrix ke can be obtained by numerically carrying out the integrals in Eq. (7.57) using the Gauss integration scheme of Eq. (7.59). 2 x 2 Gauss points shown in Figure 7.9(b) are sufficient to obtain the exact solution for the stiffness matrix given by Eq. (7.57). This is because the entry in the strain matrix, B, is a linear function of ξ or η. The integrand in Eq. (7.57) consists of BrcB, which implies multiplications of two linear functions, and hence this becomes a quadratic function. In Table 7.1, having two Gauss points sampled in each direction is sufficient to obtain the exact results for a polynomial function in that direction of order up to 3. Figure 7.9(a) and (c) show some other different and possible number of integration points in a square region.

To obtain the element mass matrix me, we substitute Eq. (7.56) into Eq. (3.75) to obtain

 

 

 

Figure 7.9. Integration points forI    in a square region.

Upon evaluation of the integral, after substitution of Eq. (7.50) into Eq. (7.60), the element mass matrix me is obtained explicitly as

In obtaining element mij in the mass matrix, the following integral has been carried out and repeatedly used:

For example, in calculating m33, we use the above equation to obtain

In practice, the integrals in Eq. (7.60) are often calculated numerically using the Gauss integration scheme.

The nodal force vector for a rectangular element can be obtained by using Eqs. (3.78), (3.79) and (3.81). Suppose the element is loaded by a distributed force fs on edge 2-3 of the element, as shown in Figure 7.8; the nodal force vector becomes

If the load is uniformly distributed within the element, and fsx and fsy are constant, the above equation becomes

where b is the half length of the side 2-3. Equation (7.65) suggests that the evenly distributed load is divided equally onto nodes 2 and 3.

The stiffness matrix ke, mass matrix me and nodal force vector fe can be used directly to assemble the global FE equation, Eq. (3.96). Coordinate transformation is needed if the orientation of the local natural coordinate does not coincide with that of the global coordinate system. In such a case, quadrilateral elements are often used, which is to be developed in the next section.

Linear Quadrilateral Elements

Though the rectangular element can be very useful, and is usually more accurate than the triangular element, it is difficult to use it for problems with any geometry rather than rectangles. Hence, its practical application is very limited. A much more practical and useful element would be the so-called quadrilateral element, that can have unparalleled edges. However, there can be a problem for the integration of the mass and stiffness matrices for a quadrilateral element, because of the irregular shape of the integration domain. The Gauss integration scheme cannot be implemented directly with quadrilateral elements. Therefore, what is required first is to map the quadrilateral element into the natural coordinates system to become a square element, so that the shape functions and the integration method used for the rectangular element can be utilized. Hence, key in the development of a quadrilateral element is the coordinate mapping. Once the mapping is established, the rest of the procedure is exactly the same as that used for formulating the rectangular element in the previous section.

Coordinate Mapping

Figure 7.10 shows a 2D domain with the shape of an airplane wing. As you can imagine, dividing such a domain into rectangular elements of parallel edges is impossible. The job can be easily accomplished by the use of quadrilateral elements with four straight but unparallel edges, as shown in Figure 7.10. In developing the quadrilateral elements, we use the same coordinate mapping that was used for the rectangular elements in the previous section. Due to the slightly increased complexity of the element shape, the mapping will become a little more involved, but the procedure is basically the same.

Consider now a quadrilateral element with four nodes numbered 1, 2, 3 and 4 in a counter-clockwise direction, as shown in Figure 7.11. The coordinates for the four nodes are indicated in Figure 7.11(a) in the physical coordinate system. The physical coordinate system can be the same as the global coordinate system for the entire structure. As there are two DOFs at a node, a linear quadrilateral element has a total of eight DOFs, like the rectangular element. A local natural coordinate system (ξ, η) with its origin at the centre of the squared element mapped from the global coordinate system is used to construct the shape functions, and the displacement is interpolated using the equation

Equation (7.66) represents a field variable interpolation using the nodal displacements. Using a similar concept, we can also interpolate the coordinates x and y themselves. In other words, we let coordinates x and y be interpolated from the nodal coordinates using the shape functions which are expressed as functions of the natural coordinates. This coordinate interpolation is mathematically expressed as

Figure 7.10. 2D domain meshed by quadrilateral elements.

Figure 7.11. Coordinates mapping between coordinate systems.

 

where X is the vector of the physical coordinates,

and N is the matrix of shape functions given by Eq. (7.50). In Eq. (7.67), xe is the physical coordinates at the nodes of the element, given by

Equation (7.67) can also be expressed explicitly as

where Ni is the shape function defined for the rectangular element in Eqs. (7.53) or (7.54). Note that, due to the unique property of the shape function, the interpolation at these nodes will be exact. For example, substituting ξ = I and η = -1 in Eq. (7.70) gives x = x2 and y = y2, as shown in Figure 7.11. Physically, this means that point 2 in the natural coordinate system is mapped to point 2 in the physical coordinate system, and vice versa. The same can be easily observed also for points 1, 3 and 4.

Let us now analyse this mapping more closely. Substituting ξ = 1 into Eq. (7.70) gives

or

Eliminating η from the above two equations gives

which represents a straight line connecting the points (x2, y2) and (x3, y3). This means that edge 2-3 in the physical coordinate system is mapped onto edge 2-3 in the natural coordinate system. The same can be observed for the other three edges. Hence, we can see that the four straight edges of the quadrilateral in the physical coordinate system correspond to the four straight edges of the square in the natural coordinate system. Therefore, the full domain of the quadrilateral element is mapped onto a square one.