<<

. 63
( 67 .)



>>

Ni = 1,
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)

<<

. 63
( 67 .)



>>