### Strain Matrix

**After mapping is performed for the coordinates,** we can now evaluate the strain matrix B. To do so in this case, it is necessary to express the differentials in terms of the natural coordinates, since the relationship between the x and y coordinates and the natural coordinates is no longer as straightforward as in the case for rectangular elements. Utilizing the chain rule in application to partial differentiation, we have

The above equations can be written in the matrix form

where J is the Jacobian matrix defined by

We now substitute the interpolation of the coordinates defined by Eq. (7.70) into the above equation, and obtain

Rewriting Eq. (7.75) to obtain

which gives the relationship between the differentials of the shape functions with respect to x and y with those with respect to ξ and η. We can now use the equation B = LN to compute the strain matrix B, by replacing all the differentials of the shape functions with respect to x and y with those with respect to ξ and η, using Eq. (7.78). This process needs to be performed numerically by a computer.

### Element Matrices

**Once the strain matrix B has been obtained,** we can proceed to evaluate the element matrices. The stiffness matrix can be obtained by Eq. (7.39). To evaluate the integration, the following formula, which has been proven by Murnaghan (1951), is used:

where det | J| is the determinate of the Jacobian matrix. Hence, the element stiffness matrix can be written as

The above integrals can then be evaluated using the Gauss integration scheme discussed in the previous section. Notice how the coordinate mapping enables us to use the Gauss integration scheme over a simple squared area.

**The shape function defined by Eq. (7.53)** is a bilinear function of ξ and η. The elements in the strain matrix B are obtained by differentiating these bilinear functions with respect to ξ and η, and by dividing by the Jacobian matrix whose elements are also bilinear functions. Therefore, the elements of BtcB det |J| are fraction functions, which cannot usually be expressed by polynomials. This means that the stiffness matrix may not be able to be evaluated exactly using the Gauss integration scheme, unlike the case for the rectangular element.

**The element mass matrix me can also be evaluated in the same way as for the rectangular element using Eq. (7.60):**

The element force vector is obtained in the same way as described for the rectangular element, because the integration for calculating the force vectors is one-dimensional line integrations. Having obtained the element matrices, the usual method of assembling the element matrices is carried out to obtain the global matrices.

### Remarks

The shape functions used to interpolate the coordinates in Eq. (7.70) are the same as those used for interpolation of the displacements. Such an element is called an isoparametric element. However, the shape functions for coordinate and displacement interpolations do not necessarily have to be the same. Using different shape functions for coordinate and displacement interpolations will lead to the development of what is known as subparametric or superparametric elements. These elements have been studied in academic research, but less often used in practical applications. Details of such elements will not be covered in this topic.

## Higher Order Elements

### Triangular Element Family General formulation of shape functions

In developing higher order elements, we make use of the area coordinate system. Figure 7.12 shows a general triangular element of order p that has nd nodes calculated by

Node i (I, J, K) is located at the I th node in the L1 direction, at the Jth node in the L2 direction, and at the Kth node in the L3 direction. From Figure 7.12, we have at any node that

The shape function can be written in the form (Argyris, 1968)

**Figure 7.12. Triangular element of order p defined under the area coordinate system. Node i (I, J, K) is located at the /th node in the L| direction, at the Jth node in the Lj direction, and at the Kth node in the L3 direction. At any node, we have I + J + K = p.**

whereare defined by Eq. (4.82), but the coordinate ξ is replaced by the area coordinates, i.e.

where a = 1, 2, 3; β = 1,J,K. For example, when a = I and β = I ,we have

Since

it is easy to confirm the delta function property of

From Eqs. (7.84) and (7.85), the order of the shape function can be found to be the same as

Since La is linear function of x and y, the order of the shape function will be

### Quadratic triangular elements

**Consider a quadratic triangular element shown** in Figure 7.13. The element has six nodes: three corner nodes and three mid-side nodes. Using Eqs (7.84) and (7.85), the shape functions can be obtained very easily. Here we demonstrate the calculation of N1. Note that for the element shown in Figure 7.13, the area coordinate L1 has three coordinate values:

**Figure 7.13. Quadratic triangular element.**

**Using Eq. (7.84),** we have

For the other two corner nodes 2, 3, we should have exactly the same equation:

For the mid-side node 4, we have

This equation is also valid for the other two mid-nodes, and therefore we have

### Cubic triangular elements

For the cubic triangular element shown in Figure 7.14 that has nine nodes, the shape function can also be obtained using Eq. (7.84), as well as four area coordinate values of (taking L1 as an example)

**We omit the process and list the results below. The reader is encouraged to confirm the results. For corner nodes (1, 2, and 3):**

**For side nodes (4-9):**

**For the interior node (10):**

**Figure 7.14. Cubic triangular element.**

**Figure 7.15. Rectangular element of arbitrary high orders.**

### Rectangular Elements Lagrange type elements

**Considering a rectangular element with**■ nodes, shown in Figure 7.15.

The element is defined in the domain ofin the natural coordinates ξ and η. Due to the regularity of the nodal distribution along both the ξ and η directions, the shape function of the element can be simply obtained by multiplying one-dimensional shape functions with respect to the ξ and η directions using the Lagrange interpolants defined in Eq. (4.82) (Zienkiewicz et al., 2000):

Due to the delta function proper of the 1D shape functions given in Eq. (4.83), it is easy to confirm that the Ni given by Eq. (7.100) is also of the delta function property.

**Figure 7.16. Nine-node rectangular element.**

**Using Eqs. (7.100) and (4.82),** the nine-node quadratic element shown in Figure 7.16 can be given by

**From Eq. (7.101),** it can easily be seen that all the shape functions are formed using the same set of nine basis functions:

which are linearly-independent.In addition, because the basis functions also contain the linear basis functions, these shape functions can also be expected to have linear field reproduction (Lemma 3), at least. Hence, they satisfy the sufficient requirements for FEM shape functions. Any other high order Lagrange type of rectangular element can be created in exactly the same way as for the nine-node element.

### Serendipity type elements

**The method used in constructing the Lagrange type of elements** is very systematic. However, the Lagrange type of elements is not very widely used, due to the presence of the interior nodes. Serendipity type elements are created by inspective construction methods. We intentionally construct high order elements without interior nodes.

**Figure 7.17. High order serendipity element. (a) Eight-node element. (b) 12-node element.**

**Consider the eight-node element shown in Figure 7.17a.** The element has four corner nodes and four mid-side nodes. The shape functions in the natural coordinates for the quadratic rectangular element are given as

whereare the natural coordinates of node j. It is very easy to observe that the shape functions possess the delta function property. The shape function is constructed by simple inspections making use of the shape function properties. For example, for the corner node 1 (wherethe shape function Ni has to pass the following three lines as shown in Figure 7.18 to ensure its vanishing at remote nodes:

**The shape N1 can** then immediately be written as

where C is a constant to be determined using the condition that it has to be unity at node 1 atwhich gives

We finally have

which is the first equation in Eq. (7.103) for j = 1.

**Figure 7.18. Construction of eight-node serendipity element. Three straight lines passing through the remote nodes of node 1 are used.**

**Figure 7.19. Construction of eight-node serendipity element. Three straight lines passing through the remote nodes of node 5 are used.**

**Shape functions at all the other corner nodes** can be constructed in exactly the same manner. As for the mid-side nodes, say node 5, we enforce the shape function passing through the following three lines as shown in Figure 7.19:

**The shape N5 can then immediately be written as**

where C is a constant to be determined using the condition that it has to be unity at node 5 atwhich gives

We finally have

which is the second equation in Eq. (7.103) for j = 5.

**Because the delta functions property is** used for the shape functions given in Eq. (7.103), they of course possess the delta function property. It can be easily seen that all the shape functions can be formed using the same set of basis functions

that are linearly-independent. From Lemmas 2 and 3, we confirm that the shape functions are partitions of unity, and at least linear field reproduction. Hence, they satisfy the sufficient requirements for FEM shape functions.

**Following a similar procedure**, the shape functions for the 12-node cubic element shown in Figure 7.17b can be written as

The reader is encouraged to figure out what lines should be used to form the shape functions listed in Eq. (7.113). When η = ηi = 1, the above equations reduce to one-dimensional cases of quadratic and cubic elements defined by Eqs. (4.84) and (4.85).