<<

. 62
( 67 .)



>>

1
N8 = (1 ’ ξ ) (1 ’ ·2 ) . (17.35)
2
Other examples may be found in Zienkiewicz [18] and Hughes [10].



17.5 Numerical integration

’ I be some function, and assume that the integral:
Let f : R
e

f ( x) dx, (17.36)
e

over the domain e of an element is to be computed. In ¬nite element computa-
tions there is a mapping from the x-space to the ξ -space, such that (see Section
16.7, on isoparametric elements)
1 dx
f ( x) dx = f ( ξ) ( ξ ) dξ . (17.37)

’1
e

φ(ξ )

This integral can be approximated with a numerical integration rule:
nint
1
g( ξ ) dξ ≈ g( ξi ) Wi , (17.38)
’1 i=1

where ξi denotes the location of an integration point and Wi the associated weight
factor.
In Fig. 17.9 an interpretation is given of the above numerical integration rule.
At a discrete number of points ξi within the interval ξ ∈[ ’1, +1] the function
value g( ξi ) is evaluated. Related to each point ξi a rectangle is de¬ned with height
g( ξi ) and width Wi . Note that it is not necessary that the point ξi is located on the
symmetry line of the rectangle. By adding up the surfaces g( ξi ) Wi of all rectangles
an approximation is obtained of the total surface underneath the function, which
is the integral. It is clear that the weight factor Wi in Eq. (17.38) can be interpreted
as the width of the interval around ξi .
The integration rule that is mostly used is the Gaussian quadrature. In that
case the locations of the integration points and weight factors are chosen so as
Shape functions and numerical integration
306

Table 17.1 Gaussian quadrature up to nint =3.
ξi
nint Wi
ξ1 = 0 W1 = 2
1
’1 1
ξ1 = √ , ξ2 = √ W1 = W2 = 1
2
3 3

5 8
3 3
ξ1 = ’ 5 , ξ2 = 0, ξ3 = W1 = W3 = 9 , W2 = 9
3 5




g(ξ)

Wi




ξi
“1 1
ξ

Figure 17.9
Numerical integration of a function φ( ξ ).


to obtain optimal accuracy for polynomial expressions of g( ξ ). In Table 17.1 the
location of the integration points and the associated weight factors are given up to
nint = 3. For two-dimensional problems the above generalizes to
1 1
= f ( x( ξ , ·) , y( ξ , ·) ) j( ξ , ·) dξ d·
f ( x, y) d
’1 ’1
e
1 1
= g( ξ , ·) dξ d·, (17.39)
’1 ’1
with j( ξ , ·) according to Eq. (16.50) and
nint nint
1 1
g( ξ , ·) dξ d· ≈ g( ξi , ·j ) Wi Wj . (17.40)
’1 ’1 i=1 j=1

The above integration scheme can be elaborated for the 9-node rectangular
Lagrangian element in Fig. 17.10 using Table 17.2.
As was already remarked in Section 17.3 integration over a triangular domain is
not trivial. For a triangular element that is formed by degeneration from a quadri-
lateral element the integration can be performed in the same way, with the same
integration points, as given above.
In the case that triangular coordinates are used, the evaluation of the integrals is
far from trivial. If the triangle coordinates »1 and »2 are maintained by eliminating
307 17.5 Numerical integration

Table 17.2 Integration points in the square and associated weight
factors.
nint point location of the integration points Wi
ξ ·
’0.7745966692 ’0.7745966692
9 1 0.3086420047
’0.7745966692
2 0.7745966692 0.3086420047
3 0.7745966692 0.7745966692 0.3086420047
’0.7745966692
4 0.7745966692 0.3086420047
’0.7745966692
5 0 0.4938271818
6 0.7745966692 0 0.4938271818
7 0 0.7745966692 0.4938271818
’0.7745966692
8 0 0.4938271818
9 0 0 0.7901234686




4 7 3
·


ξ
8 9 6

1 5 2


Figure 17.10
Position of the integration points in the square.




2
1
»2




3 1
0
0 1
»1

Figure 17.11
Mapping of triangle in a »1 , »2 -coordinate system.
Shape functions and numerical integration
308

Table 17.3 Integration points in the triangles.
nint point location of the integration points weight factors
»1 »2 »3 Wi
3 1 0.5 0.5 0 0.16667
2 0 0.5 0.5 0.16667
3 0.5 0 0.5 0.16667

7 1 0.3333333333 0.3333333333 0.3333333333 0.11250
2 0.0597158717 0.4701420641 0.4701420641 0.00662
3 0.4701420641 0.0597158717 0.4701420641 0.00662
4 0.4701420641 0.4701420641 0.0597158717 0.00662
5 0.7974269853 0.1012865073 0.1012865073 0.00630
6 0.1012865073 0.7974269853 0.1012865073 0.00630
7 0.1012865073 0.1012865073 0.7974269853 0.00630




»2
»2



6

2
2 1 4
1
7
»1 »1
3 5

3

Figure 17.12
Position of the integration points in the triangles.




the coordinate »3 these can represented in a rectangular coordinate system as given
in Fig. 17.11. It can easily be seen that the surface integral of a function φ( »1 , »2 )
can be written as
1 1’»1
φ( »1 , »2 ) d»2 d»1 . (17.41)
0 0

Despite a problem with a variable integral limit, it appeared possible to derive
numerical integration rules for triangles. In Table 17.3 the position of integration
points in triangular coordinates as well as weight factors are given for two higher-
order elements.
309 Exercises

Exercises

17.1 A MATLAB script to calculate the shape functions and the derivatives with
respect to the isoparametric coordinates of 4-noded isoparametric bi-linear
elements in some point of the elements may consist of the following code:
% Shape = program to calculate shape functions
xi= ;
eta= ;
x=[ , , , ]™;
y=[ , , , ]™;
N=[0.25*(1-xi)*(1-eta);
0.25*(1+xi)*(1-eta);
0.25*(1+xi)*(1+eta);
0.25*(1-xi)*(1+eta)]
dNdxi=[-0.25*(1-eta);
0.25*(1-eta);
0.25*(1+eta);
-0.25*(1+eta)]
dNdeta=[-0.25*(1-xi);
-0.25*(1+xi);
0.25*(1+xi);
0.25*(1-xi)]


y


1 (10, 6)
(5, 6)

32
2 4
(0, 4)
(8, 5)
4
1
23
(6, 2)
1

(0, 0)
x



To study some of the properties of the shape functions the script will be
completed by considering the elements in the ¬gure above and extended
in the following. In this ¬gure the local node numbering has also been
indicated.
Shape functions and numerical integration
310

(a) Use the program to test whether
n

<<

. 62
( 67 .)



>>