i=1

and

n n

‚Ni ‚Ni

= = 0,

‚ξ ‚·

i=1 i=1

for a number of combinations of ( ξ , ·).

(b) Extend the program to calculate the Jacobian matrix:

‚x ‚x

‚ξ ‚·

x,ξ = ,

‚y ‚y

‚ξ ‚·

and the Jacobian determinant j = det( x,ξ ).

(c) Determine the Jacobian determinant in the following points:

( ξ1 , ·1 ) = (0, 0)

( ξ2 , ·2 ) = (0.5, 0.5)

( ξ3 , ·3 ) = (1, 0)

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

for both elements which are shown in the ¬gure. What can you con-

clude from this?

17.2 Consider the element that is given in the ¬gure below. The element shape

functions of this element are derived after degeneration of a 4-noded

quadrilateral element by coalescence of two nodes in the same way as

discussed in Section 17.3.

y

(0, 1) 3

2

1

x

(x1, 0)

(0, 0)

Give analytical expressions for the shape functions Ni ( ξ , ·).

(a)

Compute the Jacobian determinant for ξ = · = 0.

(b)

Plot the result as a function of x1 on the interval [ ’2, +2] and

(c)

comment on the result.

311 Exercises

17.3 Consider the one-dimensional 4-noded element as given in the ¬gure.

ξ = “1 ξ = “1/3 ξ = +1/3 ξ = +1

Use Lagrange polynomials to derive the element shape functions.

17.4 After solving the stationary convection-diffusion equation for the 4-noded

quadrilateral element e the following element solution vector has been

found:

¤

⎡

0.4

⎢ 0.1 ⎥

⎥

⎢

ue = ⎢ ⎥.

⎣ 0.3 ¦

∼

1

(a) Use the script from Exercise 17.1 as the basis to calculate the uh at

the point ( ξ , ·) = ( 0.5, 0.5).

(b) The nodal coordinates of the element are given:

⎡ ¤

5 2

⎢ ⎥

⎢ ⎥

15 1

coord = ⎢ ⎥.

⎣ ¦

10 4

4 4

Calculate ‚uh /‚x and ‚uh /‚y for ( ξ , ·) = ( 0.5, 0.5).

17.5 Consider the quadratic 6-node triangular element in the ¬gure.

y

3

7

6

5

5

6

4

2

3

4

2

1

1

x

3 5 6

2 7

1 4

The solution vector for this element after a computation is given as

Shape functions and numerical integration

312

⎡ ¤

2

⎢ ⎥

⎢ ⎥

5

⎢ ⎥

⎢ ⎥

3

ue = ⎢ ⎥.

⎢ ⎥

∼

2

⎢ ⎥

⎣ ¦

0

2

Write a MATLAB program to calculate the solution uh at point

(a)

( x, y) = ( 4, 4).

Give expressions for ‚uh /‚x and ‚uh /‚y as a function of x and y.

(b)

18 Infinitesimal strain elasticity

problems

18.1 Introduction

One of the ¬rst applications of the Finite Element Method in biomechanics has

been the analysis of the mechanical behaviour of bone [2]. In particular the impact

of prosthesis implants has been investigated extensively. An example of a ¬nite

element mesh used to analyse the mechanical stresses and strain in a human femur

is given in Fig. 18.1(a). In Fig. 18.1(b) the femur head is replaced by a prosthesis.

In this case the prosthesis has different mechanical properties compared to the

bone, leading to high stress concentrations at some points in the bone and stress

shielding (lower stresses than normal) in other parts. This normally leads to a

remodelling process in the bone that has to be accounted for when new prostheses

are designed. The purpose of this chapter is to introduce the ¬nite element theory

that forms the basis of these analyses.

18.2 Linear elasticity

Neglecting inertia, the momentum equation, Eq. (11.9), may be written as

∇ · σ + f = 0, (18.1)

where ∇ denotes the gradient operator, σ the Cauchy stress tensor and f = ρq a

given distributed volume load. This equation should hold at each position within

the domain of interest , having boundary , and must be supplemented with

suitable boundary conditions. Either the displacement ¬eld u is speci¬ed along u :

u = u0 at u, (18.2)

or the external load along is speci¬ed:

p

σ · n = p at p. (18.3)

The vector n denotes the unit outward normal at the boundary .

Infinitesimal strain elasticity problems

314

Hip implant

z z

y y

x x

(a) (b)

Figure 18.1

(a) Mesh to analyse stress and strain distributions in a human femur (b) a ¬nite element mesh of a

femur, where the femur head has been replaced by a prosthesis. (Image courtesy of Mr Bert van

Rietbergen.)

If isotropic, linearly elastic material behaviour according to Hooke™s law

is assumed, see Section 12.2, and the Cauchy stress tensor is related to the

in¬nitesimal strain tensor µ by

σ = Ktr(µ) I + 2Gµ d , (18.4)

where K is the compression modulus and G the shear modulus. The in¬nitesimal

strain tensor is given by:

1

µ= ( ∇u + ( ∇u)T ) , (18.5)

2

where u denotes the displacement ¬eld. Furthermore, tr( µ) is the trace and µd

the deviatoric part of the in¬nitesimal strain tensor. The shear modulus G and the

compression modulus K may be expressed in terms of Young™s modulus E and

Poisson™s ratio ν via

E

G= (18.6)

2(1 + ν)

E

K= . (18.7)

3(1 ’ 2ν)

Introduce the Cartesian basis {ex , ey , ez }. With respect to this basis the Cauchy

stress tensor may be written as

σ = σxx ex ex + σxy ex ey + · · · + σzz ez ez , (18.8)

315 18.3 Weak formulation

or, using the Einstein summation convention (the presence of double indices

implies summation over these indices):

σ = σij ei ej , (18.9)

with i, j = x, y, z. The Cauchy stress matrix is given by

⎡ ¤

σxx σxy σxz

⎢ ⎥

σ = ⎣ σxy σyy σyz ¦ , (18.10)

σxz σyz σzz

where use has been made of the symmetry of the Cauchy stress tensor σ T = σ .

The strain matrix is given by

⎡ ¤

µxx µxy µxz

⎢ ⎥

µ = ⎣ µxy µyy µyz ¦ . (18.11)

µxz µyz µzz

The trace of the strain tensor is given by

tr( µ) = µxx + µyy + µzz , (18.12)

and the deviatoric part of the strain matrix is given by

⎡ ¤

3 µxx ’ 3 ( µyy + µzz ) µxy µxz

2 1

⎢ ⎥

µ =⎣

d

µxy 3 µyy ’ 3 ( µxx + µzz ) µyz

2 1

¦.

µxz µyz 3 µzz ’ 3 ( µxx + µyy )

2 1

(18.13)

18.3 Weak formulation

As before, the weak formulation based on the method of weighted residuals is

obtained by pre-multiplying Eq. (18.1) with a weight function. Since Eq. (18.1)