2.3 Numerical Integration

Back

To obtain the force, the stiffness matrix and the mass matrix (see Equations (2.22)–(2.24)),

integration over the element is required. To calculate an integral of the form

I = _

V0e

f (X) dV (2.104)

numerical integration is used. Indeed, f (X) is frequently a complicated function of space (it

usually contains material properties depending on the temperature, which is itself a function

of space) and the shape of an element can be quite irregular. Consequently, analytical

integration is not feasible. Before applying numerical integration to Equation (2.104) the

integration domain is transformed from the global to the local element coordinate system:

I = _

V0eL

f [X(ξ, η, ζ)]J

∗

(ξ, η, ζ) dξ dη dζ =: _

V0eL

g(ξ, η, ζ ) dξ dη dζ (2.105)

where J

∗

(ξ, η, ζ) is the Jacobian determinant of the transformation X(γ ) (see

Equation (2.17)). In this way, the integration domain is identical for all elements belonging

to the same type, for example, brick elements. Now, the analytical integration is approximated

by a numerical integration scheme (Stroud 1971). This basically means that the

integral in Equation (2.105) is replaced by a linear combination of function values at

specific locations, the so-called integration points, within the domain of integration:

_

V0eL

g(ξ, η, ζ ) dξ dη dζ ≈

N

_

i=1

g(ξi, ηi, ζi)wi . (2.106)

76 LINEAR MECHANICAL APPLICATIONS

The value of the weights, the location of the integration points and their number constitute

together an integration scheme. Different schemes lead to different calculational expenditure

and different accuracy. For finite element calculations, the Gauss schemes are very popular,

because of their high accuracy compared to the numerical expenditure. Since the integration

schemes essentially depend on the shape of the domain, a distinction is made between

hexahedral, tetrahedral and wedge elements.

2.3.1 Hexahedral elements

The domain in local coordinates for a hexahedral element is a cube extending from −1 to +1

(−1 ≤ ξ, η, ζ ≤ 1) along each coordinate axis. The integration schemes are symmetric in

each direction. The lowest scheme has one integration point in each direction (1 × 1 × 1 =

1, Figure 2.7), the next ones have two (2 × 2 × 2 = 8, Figure 2.8) or three (3 × 3 × 3 = 27,

Figure 2.9) points in each direction. The location of the integration points and their weights

are summarized in Table 2.1.

The 1 × 1 × 1 scheme, the 2 × 2 × 2 scheme and the 3 × 3 × 3 scheme are exact for

a constant function, a trilinear function and a triquadratic function respectively. Therefore,

the 2 × 2 × 2 scheme represents full integration for a linear element (8–node brick) and

the 3 × 3 × 3 scheme stands for full integration in a quadratic element (20–node brick).

The term reduced integration is used if one selects the next coarser scheme: 1 × 1 × 1 for

linear elements and 2 × 2 × 2 for quadratic elements. Reduced integration frequently has

a beneficial effect: it produces less shear locking and less volumetric locking; therefore, it

is ideal for plates, shells and incompressible materials (rubber, plasticity in metals). Furthermore,

Barlow has shown ((Barlow 1976), see also (Mackinnon and Carey 1989) and

1

1

2

3

4

5

6

7

8

Figure 2.7 Hexahedral element: 1 × 1 × 1 scheme

LINEAR MECHANICAL APPLICATIONS 77

1

2

3

4

5

6

7

8

1

2

3

4

5

6

7

8

Figure 2.8 Hexahedral element: 2 × 2 × 2 scheme

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

1

2

3

4

5

6

7

8

Figure 2.9 Hexahedral element: 3 × 3 × 3 scheme

78 LINEAR MECHANICAL APPLICATIONS

Table 2.1 Location of the integration points in hexahedral

elements.

Scheme Location (ξi, ηi , ζi ) Number Weight

and any perturbation

1 × 1 × 1 (0, 0, 0) 1 8

2 × 2 × 2 _± 1 √

3

,± 1 √

3

,± 1 √

3

_ 8 1

3 × 3 × 3

_

±_3

5

,±_3

5

,±_3

5

_

8 _5

9

_3

_

0,±_3

5

,±_3

5

_

12 _8

9

__5

9

_2

_

0, 0,±_3

5

_

6 _8

9

_2 _5

9

_

(0, 0, 0) 1 _8

9

_3

(Liew and Rajendran 2002)) that the reduced integration points are the so-called superconvergent

points, in which the stress is one order more accurate than in any other point.

However, because of reduced integration, so-called zero-energy modes can arise, leading

to hourglassing. Shear locking, volumetric locking and hourglassing are discussed in

Section 2.5.

2.3.2 Tetrahedral elements

For tetrahedral elements, the integration domain in local coordinates is depicted in Figure

2.10. The most frequently used Gauss integration schemes are summarized in Table 2.2.

Linear tetrahedral elements are usually integrated with one integration point, quadratic

elements with four. Figure 2.10 visualizes the scheme with four integration points. The

scheme with 15 integration points improves the condition of the consistent mass matrix (cf

Section 2.11.6). For tetrahedral elements, the term reduced integration is not used.

2.3.3 Wedge elements

The integration domain for a wedge element consists of a prism with triangular lower and

upper surface (Figure 2.11). Linear wedge elements are usually integrated with a 2-point

scheme, quadratic wedges with a 9-point scheme. The 18-point scheme is used for the

integration of the consistent mass matrix (cf Section 2.11.6). The integration schemes are

summarized in Table 2.3.

LINEAR MECHANICAL APPLICATIONS 79

1

2

3

4

1

2

3

4

Figure 2.10 Tetrahedral element: 4 integration points

Table 2.2 Location of the integration points in tetrahedral elements.

Total number of Location (ξi, ηi , ζi, 1 − ξi

− ηi

− ζi ) Number Weight

integration points and any perturbation

1 _1

4

,

1

4

,

1

4

,

1

4

_ 1

1

6

4

_5 −

√

5

20

,

5 −

√

5

20

,

5 −

√

5

20

,

5 + 3

√

5

20

_

4

1

24

15 _1

4

,

1

4

,

1

4

,

1

4

_ 1

16

810

_7 −

√

15

34

,

7 −

√

15

34

,

7 −

√

15

34

,

13 + 3

√

15

34

_

4

2665 + 14

√

15

226 800

_7 +

√

15

34

,

7 +

√

15

34

,

7 +

√

15

34

,

13 − 3

√

15

34

_

4

2665 − 14

√

15

226 800

_10 − 2

√

15

40

,

10 − 2

√

15

40

,

10 + 2

√

15

40

,

10 + 2

√

15

40

_

6

20

2268

80 LINEAR MECHANICAL APPLICATIONS

1

2

3

4

5

6

7

8

9

1

2

3

4

5

6

Figure 2.11 Wedge element: 9 integration points

Table 2.3 Location of the integration points in wedge elements.

Total number of Location (ξi, ηi, 1 − ξi

− ηi

− ζi , ζi ) Number Weight

integration points and any perturbation

2 _1

3

,

1

3

,

1

3

;± 1 √

3

_ 2

1

2

9

_1

6

,

1

6

,

4

6

;±_3

5

_

6

5

54

_1

6

,

1

6

,

4

6

; 0_ 3

8

54

18

_1

6

,

1

6

,

4

6

;±_3

5

_

6

1

12

_1

6

,

1

6

,

4

6

; 0_ 3

2

15

_1

2

,

1

2

, 0;±_3

5

_

6

1

108

_1

2

,

1

2

, 0; 0_ 3

2

135

LINEAR MECHANICAL APPLICATIONS 81

X

ξ

η

−1

−1

1

dξ 1

dη

dξ = 0

dη = 0

dX

dξ dξ

dX

dη dη

XL

Figure 2.12 Mapping a surface element in local coordinates onto global coordinates

2.3.4 Integration over a surface in three-dimensional space

Occasionally, the domain of an integral is a surface. For instance, for a distributed pressure,

the first term in Equation (2.23) takes the form

−_

A0e

pNKϕi dAe = −_

A0e

pϕi dAK

e (2.107)

or, generically,

I = _

A0e

f (X) dA. (2.108)

On the left-hand side of Figure 2.12, the surface is shown in local coordinates, on the

right-hand side in global coordinates. The infinitesimal surface dA satisfies

dA = ∂X

∂ξ

dξ × ∂X

∂η

dη (2.109)

where “×” is the vector product. Consequently, Equation (2.108) can be replaced by

I = _

AeL

f [X(ξ, η)]J

∗ dξ dη (2.110)

where

J

∗ := ∂X

∂ξ

× ∂X

∂η

= ∂XK

∂ξ

∂XL

∂η

GK × GL (2.111)

can be considered as a Jacobian vector. Since

GK × GL = eKLMGM_detG_ (2.112)

82 LINEAR MECHANICAL APPLICATIONS

Equation (2.111) can also be written as

J

∗ = eKLM

∂XK

∂ξ

∂XL

∂η

GM_detG_ (2.113)

or

J

∗ =

__________

G1 G2 G3

∂X1

∂ξ

∂X2

∂ξ

∂X3

∂ξ

∂X1

∂η

∂X2

∂η

∂X3

∂η

__________

_detG_ (2.114)

where the vertical lines denote the determinant. Notice that the vector product of two vectors

yields a one-form. This confirms the one-form nature of a differential surface element.

Equation (2.107) is the expression for the force in direction K in local node i due to

a distributed pressure on A0e. Focusing on a hexahedral element type and assuming that

x = ξ , y = η, z = 0 and that the pressure is constant, the force Fzi in z-direction in local

node i takes the form

Fzi = −p _ 1

−1

_ 1

−1

ϕi(ξ, η) dξ dη. (2.115)

The total force Fz on the surface amounts to

Fz = −4p. (2.116)

Accordingly, the relative force in local node i satisfies

Fzi/Fz = 1

4

_ 1

−1

_ 1

−1

ϕi(ξ, η) dξ dη. (2.117)

The shape functions in the face are the three-dimensional shape functions from

Section 2.2 for which one local coordinate is kept constant. Performing the integration

in Equation (2.117) and similarly for the other element types leads to the force distributions

in Figure 2.13. One notices that for linear elements each node takes the same amount

of force. This is not the case for quadratic elements. For 10-node tetrahedral elements, the

vertex nodes do not take any force at all and the middle nodes carry the complete force.

For 20-node brick elements, the middle nodes carry even more than the complete force

resulting in tensile forces in the vertex nodes. This substantially complicates the detection

of contact conditions in quadratic elements.

Навигация