x x

x

∼3 ∼3

∼3

x x x

∼6 x

∼5 ∼6 ∼5

x

∼7

x

x x

x

x x x

x ∼2

∼1 ∼4

∼1 ∼1 ∼4 ∼2

∼2

a. 3 nodes b. 6 nodes c. 7 nodes

Figure 17.4

Triangular elements with 3, 6 or 7 nodes.

»3

y

3

3

»2

2

1

2 1

x

»1

Figure 17.5

The mapping of a triangle from the global coordinate system to the local coordinate system with

triangle coordinates.

301 17.3 Linear triangular element

(a) 3-node triangular element, linear interpolation:

N1 = »1

N2 = »2

N3 = »3 (17.21)

(b) 6-node triangular element, quadratic interpolation:

N1 = »1 (2»1 ’ 1)

N2 = »2 (2»2 ’ 1)

N3 = »3 (2»3 ’ 1)

N4 = 4»1 »2

N5 = 4»2 »3

N6 = 4»1 »3 (17.22)

(c) 7-node triangular element, bi-quadratic interpolation:

N1 = »1 (2»1 ’ 1) + 3»1 »2 »3

N2 = »2 (2»2 ’ 1) + 3»1 »2 »3

N3 = »3 (2»3 ’ 1) + 3»1 »2 »3

N4 = 4»1 »2 ’ 12»1 »2 »3

N5 = 4»2 »3 ’ 12»1 »2 »3

N6 = 4»1 »3 ’ 12»1 »2 »3

N7 = 27»1 »2 »3 . (17.23)

The factor »1 »2 »3 is called a ˜bubble™ function giving zero contributions along the

boundaries of the element.

There are two major differences between the method used in Section 17.2 and

the method with triangle coordinates:

• Determining the derivatives of these shape functions with respect to the global coor-

dinates is not trivial. Consider the derivative of a shape function Ni with respect to x.

Applying the chain rule for differentiation would lead to

‚Ni ‚Ni ‚»1 ‚Ni ‚»2 ‚Ni ‚»3

= + + . (17.24)

‚x ‚»1 ‚x ‚»2 ‚x ‚»3 ‚x

But, by de¬nition a partial derivative as to one variable implies that the other variables

have to be considered as constant. In this case a partial derivative with respect to »1

means that, when this derivative is determined, »2 and »3 have to be considered constant.

However, the »i ™s are related by Eq. (17.17). Only two variables can be considered as

Shape functions and numerical integration

302

independent. A solution for this dilemma is to eliminate »3 from the shape functions by

using:

»3 = 1 ’ »1 ’ »2 . (17.25)

A more obvious way of solving the problem is to substitute Eq. (17.18) into the equa-

tions for the shape functions, thus eliminating all triangular coordinates and directly

determine ‚Ni /‚x and ‚Ni /‚y.

• The second issue is the difference in integration limits which have to correspond with

a triangle. For the square element in Section 17.2 the domain for integration is simple,

meaning that the surface integral can be split into two successive single integrals with

ξ and · as variables. The integration limits are ’1 to +1. For the triangle this is more

complicated and the limits of integration now involve the coordinate itself. This item

will be discussed shortly in Section 17.5.

17.4 Lagrangian and Serendipity elements

In principle higher-order elements are more accurate than the linear ones discussed

so far. However, the computation of element coef¬cient matrices and element

arrays is more expensive for higher-order elements, and the cost-effectiveness

depends on the particular problem investigated. Cost-effectiveness in the sense

that there is a trade-off between the accuracy, using a smaller number of

higher-order elements, versus using more linear elements.

Higher-order elements can systematically be derived from Lagrange polynomi-

als. A (one-dimensional) set of Lagrange polynomials on an element with domain

ξ1 ¤ ξ ¤ ξn is de¬ned by

n

b=1,b=a ( ξ ’ ξb )

la ( ξ ) =

n’1

(17.26)

n

b=1,b=a ( ξa ’ ξb )

( ξ ’ ξ1 ) . . . ( ξ ’ ξa’1 ) ( ξ ’ ξa+1 ) . . . ( ξ ’ ξn )

= ,

( ξa ’ ξ1 ) . . . ( ξa ’ ξa’1 ) ( ξa ’ ξa+1 ) . . . ( ξa ’ ξn )

with n the number of nodes of the element and with a = 1, 2, . . . n referring to a

node number. Notice that the above polynomial is of the order ( n ’ 1).

For instance, ¬rst-order (linear) polynomials are found for n = 2, hence

ξ ’ ξ2

l1 =

1

, (17.27)

ξ1 ’ ξ2

and

ξ ’ ξ1

l2 =

1

. (17.28)

ξ2 ’ ξ1

303 17.4 Lagrangian and Serendipity elements

17.4.1 Lagrangian elements

The shape functions of an element of order ( n ’ 1) in one dimension are chosen as

Na = la .

n’1

(17.29)

The quadratic shape function associated with node 1 of the element depicted in

Fig. 17.6 (with ξ1 = ’1, ξ2 = 0, ξ3 = 1) satis¬es

( ξ ’ ξ2 ) ( ξ ’ ξ3 ) 1

N 1 ( ξ ) = l1 ( ξ ) = = ξ ( ξ ’ 1) .

2

(17.30)

( ξ1 ’ ξ2 ) ( ξ1 ’ ξ3 ) 2

Likewise it follows that

N2 ( ξ ) = ’( ξ + 1) ( ξ ’ 1) = 1 ’ ξ 2 , (17.31)

and

1

N3 ( ξ ) =

ξ ( ξ + 1) . (17.32)

2

In two dimensions, the shape functions for the 9-node element as visualized

in Fig. 17.7 are formed by multiplication of two Lagrangian polynomials.

Leading to:

1

N1 ( ξ , ·) = ξ ( ξ ’ 1) ·( · ’ 1)

4

1

N2 ( ξ , ·) = ’ ( ξ + 1) ( ξ ’ 1) ·( · ’ 1)

2

ξ

1 2 3

Figure 17.6

A one-dimensional quadratic element, ’1 ¤ ξ ¤ 1.

·

7 6 5

8 9 4 ξ

1 2 3

Figure 17.7

A two-dimensional quadratic element, ’1 ¤ ξ ¤ 1, ’1 ¤ · ¤ 1.

Shape functions and numerical integration

304

1

N3 ( ξ , ·) = ’ ξ ( ξ + 1) ·( · ’ 1)

4

1

N4 ( ξ , ·) = ’ ξ ( ξ + 1) ( · + 1) ( · ’ 1)

2

1

N5 ( ξ , ·) = ξ ( ξ + 1) ( · + 1) ·

4

1

N6 ( ξ , ·) = ’ ( ξ + 1) ( ξ ’ 1) ( · + 1) ·

2

1

N7 ( ξ , ·) = ’ ξ ( ξ ’ 1) ( · + 1) ·

4

1

N8 ( ξ , ·) = ’ ξ ( ξ ’ 1) ( · + 1) ( · ’ 1)

2

N9 ( ξ , ·) = ( ξ + 1) ( ξ ’ 1) ( · + 1) ( · ’ 1) . (17.33)

17.4.2 Serendipity elements

For serendipity elements no internal nodes are used. Consider a ˜quadratic ele-

ment™, as depicted in Fig. 17.8 (right). The shape functions of the corner nodes are

de¬ned by

1

N1 =(1 ’ ξ ) (1 ’ ·) ( ’ ξ ’ · ’ 1)

4

1

N2 = (1 + ξ ) (1 ’ ·) ( + ξ ’ · ’ 1)

4

1

N3 = (1 + ξ ) (1 + ·) ( + ξ + · ’ 1)

4

1

N4 = (1 ’ ξ ) (1 + ·) ( ’ ξ + · ’ 1) , (17.34)

4

while the shape functions for the mid-side nodes read

Lagrange Serendipity

7 6 5 4 7 3

8 9 4 8 6

1 5 2

1 2 3

Figure 17.8

Example of a Serendipity element compared to the ˜equal order™ Lagrangian element.

305 17.5 Numerical integration

1

N5 =

(1 ’ ξ 2 ) (1 ’ ·)

2

1

N6 = (1 + ξ ) (1 ’ ·2 )

2

1

N7 = (1 ’ ξ 2 ) (1 + ·)

2