. 46
( 67 .)


σyy = 0
σxy = ’ y2 ,

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

σ = ’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

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



z x


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


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

+ 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

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:

= P at
c p, (14.3)
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

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

glucose t = 12 [h] 2.5

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
+ f = 0,
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
w( x) g( x) dx = 0 for all w, (14.5)

and to emphasize this important equivalence:
g( x) = 0 a¤x¤b” w( x) g( x) dx = 0
on for all w. (14.6)

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-


. 46
( 67 .)