b2

σxy = ’ y2 ,

4

with c1 and c2 constants.

Determine c2 based on the relationship between σxy and P and subse-

quently, determine c1 by means of the local equilibrium equations.

13.5 In the environment of the origin of a Cartesian xyz-coordinate system it is

given that a (two-dimensional) stationary velocity ¬eld in an incompress-

ible Newtonian ¬‚uid (with density ρ) can be described as

v = ±( ’yex + xey ) with ± a constant.

Based on this velocity ¬eld the deformation rate tensor D can be calcu-

lated: D = 0 with 0 the zero tensor. With substitution into the constitutive

equation:

σ = ’pI + 2·D,

with · the viscosity and σ the stress tensor (and with p = 0 originated by

the applied boundary conditions) it follows that σ = 0.

Then, by means of the Navier“Stokes equation it is found that the dis-

tributed load q (force per unit mass) necessary to realize the described ¬‚ow

¬eld is given by

r = xex + yey .

q = ’± 2 r with

Solution strategies for solid and fluid mechanics problems

230

Interpret the stated results: why D = 0 is obtained for the given ¬‚ow ¬eld

and what is the physical meaning of q = ’± 2 r?

13.6 An incompressible Newtonian ¬‚uid (density ρ, viscosity ·) is located in a

cavity (rectangular) with a cross section as given in the ¬gure: ’a ¤ x ¤ a,

’b ¤ y ¤ b. Because the top wall closing the cavity is translating in the

V

y

2b

z x

2a

x-direction with velocity V, a stationary (two-dimensional) ¬‚uid ¬‚ow (in

the shape of a vortex) will develop in the cavity: vx = vx ( x, y) , vy =

vy ( x, y) , vz = 0. No slip conditions hold at all walls of the cavity (moving

as well as stationary).

Which inconsistency (˜something leading to contradiction™) can be

observed in the above given problem de¬nition?

13.7 A compressible medium (for example a gas) is forced in the x-direction

through a contraction (see ¬gure). With respect to the components of the

x

velocity ¬eld v of the medium, de¬ned in the Cartesian xyz-coordinate sys-

∼

tem, it is assumed that vx only depends on x, thus vx = vx ( x), completed

with vy = vz = 0.

Why is it impossible that this assumption is correct?

13.8 In the neighbourhood of the origin of a Cartesian xyz-coordinate system

the stationary (two-dimensional) ¬‚ow ¬eld in an incompressible Newtonian

¬‚uid (with density ρ) can be described as

v = ±( ’yex + xey ) with ± a constant.

231 Exercises

With respect to the constitutive equation:

σ = ’pI + 2·D,

with · the viscosity, σ the stress tensor and D the deformation rate tensor, it

can be stated (based on the prescribed boundary conditions) that p satis¬es

p = 0.

Determine by means of the Navier“Stokes equation the distributed load q

(per unit mass) necessary to realize the given ¬‚ow pattern.

14 Solution of the one-dimensional

diffusion equation by means of

the Finite Element Method

In the present and following chapters extensive use will be made of a simple

¬nite element code mlfem_nac. This code, including a manual, can be freely

downloaded from the website: www.mate.tue.nl/biomechanicsbook.

The code is written in the program environment MATLAB. To be able to use

this environment a licence for MATLAB has to be obtained. For information about

MATLAB see: www.mathworks.com.

14.1 Introduction

It will be clear from the previous chapters that many problems in biomechanics

are described by (sets of) partial differential equations and for most problems it

is dif¬cult or impossible to derive closed form (analytical) solutions. However, by

means of computers, approximate solutions can be determined for a very large

range of complex problems, which is one of the reasons why biomechanics as

a discipline has grown so fast in the last three decades. These computer-aided

solutions are called numerical solutions, as opposed to analytical or closed form

solutions of equations. The present and following chapters are devoted to the

numerical solution of partial differential equations, for which several methods

exist. The most important ones are the Finite Difference Method and the Finite

Element Method. The latter is especially suitable for partial differential equations

on domains with complicated geometries, material properties and boundary con-

ditions (which is nearly always the case in biomechanics). That is why the next

chapters focus on the Finite Element Method. The basic concepts of the method

are explained in the present chapter. The one-dimensional diffusion equation will

be used as an example to illustrate the key features of the ¬nite element method.

Clearly, the one-dimensional diffusion problem can be solved analytically for a

wide variety of parameter choices, but the structure of the differential equation is

representative of a much larger class of problems to be discussed later. In addition,

the diffusion equation and the more extended, instationary (convection) diffusion

equation play an important role in many processes in biomechanics.

233 14.2 The diffusion equation

14.2 The diffusion equation

The differential equation that describes the one-dimensional diffusion problem is

given by

du

d

+ f = 0,

c (14.1)

dx dx

where u( x) is the unknown function, c( x) > 0 a given material characteristic func-

tion and f ( x) a given source term. This differential equation is de¬ned on a one-

dimensional domain that spans the x-axis between x = a and x = b while the

boundary, which is formally denoted by , is located at x = a and x = b.

Eq. (14.1) is an adapted form of the diffusion equation Eq. (13.43), introduced

in the previous chapter. Different symbols for the unknown (u instead of ρ) and

coef¬cient (c instead of D) are used to emphasize the general character of the

equation, applicable to different kinds of problems (see below). Furthermore, the

coef¬cient c can be a function of x and a source term f ( x) is introduced.

Two types of boundary conditions can be discerned. Firstly, the essential bound-

ary condition, which must be speci¬ed in terms of u. For the derivations that

follow the boundary at x = a is chosen, to specify this type of boundary

condition:

u = U at u, (14.2)

where u denotes the boundary of the domain at x = a. Secondly, a natural

boundary condition may be speci¬ed. Here the boundary at x = b is chosen to

specify the ¬‚ux c du/dx:

du

= P at

c p, (14.3)

dx

where p denotes the boundary at x = b. For diffusion problems, an essential

boundary condition must be speci¬ed to have a well-posed boundary value prob-

lem. This is not necessarily the case for natural boundary conditions; they may be

absent.

Example 14.1 The diffusion equation describes a large range of problems in biomechanics and

is applicable in many different areas. In the way it was introduced in Section 13.4

the unknown u represents the concentration of some matter, for example: oxy-

gen in blood, proteins or other molecules in an extracellular matrix or inside a

cell, medication in blood or tissue. In that case the term f can be either a source

term (where matter is produced) or a sink term (where matter is consumed).

Numerical solution of one-dimensional diffusion equation

234

[mMol]

glucose t = 12 [h] 2.5

2

1.5

1

construct

porous ¬lter

(a) (b)

Figure 14.1

(a) Schematic view of a bioreactor system designed for growing articular cartilage tissue (b) result

of a numerical calculation of the concentration of glucose in the tissue engineered construct. In the

analysis it is assumed that the glucose concentration in the medium surrounding the construct is

constant (essential boundary condition) and that glucose is consumed by the cells in the construct

(sink term). Because of symmetry, only the right half of the construct is modelled. Adapted

from [17].

In the human body diffusion is very important, but also in in vitro experimen-

tal set-ups in the laboratory, for example in tissue engineering applications (see

Fig. 14.1).

Example 14.2 Another completely different process which is also relevant in biomechanics is

steady state one-dimensional heat conduction with a source term:

d dT

» + f = 0,

dx dx

where » is the Fourier coef¬cient of heat conduction and f a heat source term.

Further, T( x) represents the temperature. The boundary conditions might be for-

mulated as: prescribed temperature at u and prescribed heat ¬‚ux at p . It should

be noticed that the mathematical structure of this heat conduction problem is fully

equivalent to the structure of the diffusion problem.

Example 14.3 A third example of a completely different diffusion type problem is the uniax-

ial tension or compression of a bar as introduced in Chapter 6, governed by the

equation:

du

d

+ f = 0,

EA

dx dx

where E is the Young™s modulus, A the cross section of the bar and f a dis-

tributed force per unit length. The unknown u( x) is the displacement ¬eld of the

bar (also see Eq. (6.20)). Boundary conditions may be imposed as a prescribed

displacement on u and a prescribed force at p .

235 14.3 Method of weighted residuals and weak form

The approximate solution of Eq. (14.1) is found by transforming the differential

equation into a discrete set of ordinary equations:

d du

+ f = 0 ’’ K u = f .

c (14.4)

∼ ∼

dx dx

The array u contains approximations of the continuous solution u of the differen-

∼

tial equation at a ¬nite number of locations on the x-axis. Increasing the number

of points de¬ning u should lead to an increased accuracy of the approximation of

∼

u. A particularly attractive feature of the Finite Element Method is that the spatial

distribution of these points does not need to be equidistant and can be chosen such

that accurate solutions can be obtained with a limited number of points, even on

complicated geometries (in the multi-dimensional case) or problems with large

gradients in the solution.

The ¬nite element method proceeds along three well-de¬ned steps.

(i) Transformation of the original differential equation into an integral equation by

means of the principle of weighted residuals.

(ii) Discretization of the solution u by interpolation. If an approximation of the solu-

tion u is known at a ¬nite number of points (nodes) an approximation ¬eld may be

constructed by interpolation between these point (nodal) values.

(iii) Using the discretization the integral equation is transformed into a linear set of

equations from which the nodal values u can be solved.

∼

14.3 Method of weighted residuals and weak form of the model problem

First of all the differential equation is transformed into an integral equation by

means of the weighted residuals method. Suppose that a given function g( x) = 0

on a certain domain a ¤ x ¤ b, then this formulation is equivalent to requiring

b

w( x) g( x) dx = 0 for all w, (14.5)

a

and to emphasize this important equivalence:

b

g( x) = 0 a¤x¤b” w( x) g( x) dx = 0

on for all w. (14.6)

a

The function w( x) is called the weighting function, and is assumed to be a con-

tinuous function on the integration domain. The equivalence originates from the

requirement that the integral equation must hold for all possible weighting func-