. 56
( 67 .)


u( x, t = 0) = 1 for 0.25 ¤ x ¤ 0.5.
Solve this problem in the same way.
15.4 Many biological materials can be described as a mixture of a porous solid
and a ¬‚uid, like articular cartilage [13], skin [14], intervertebral disk [6],
heart tissue [11], etc. A con¬ned compression test can be used to determine
the material parameters. A circular specimen is placed in a con¬ning ring.
A porous ¬lter supports the solid phase, while ¬‚uid from the specimen can
be expelled. On top a tight-¬tting indenter is placed to mechanically load
the specimen with a constant force. Because of the con¬ning ring the speci-
men can only deform in one direction. When a step load is applied to the

loading shaft

saline solution
tissue sample
confining ring


tissue, at ¬rst a pressure will be built up in the ¬‚uid. Immediately after load-
ing the ¬‚uid will start moving through the ¬lter and the solid will gradually
take over the load. After some time the pressure in the ¬‚uid is zero and all
load is taken by the solid. It can be derived from the theory of mixtures
The one-dimensional convection-diffusion equation

that the pressure in the ¬‚uid is described by the instationary diffusion
‚p ‚ 2p
= HK 2 ,
‚t ‚x
with K the con¬ned compression modulus (also called aggregate modulus)
and H the permeability of the porous solid. The boundary conditions are:
• p = 0 at x = 0; free drainage at the porous ¬lter
• ‚p/‚x = 0 at x = 10’2 [m]; no ¬‚ow through the surface where the
indenter makes contact with the specimen
The initial condition is p = 1000 [Nm’2 ] at every point in the specimen
at t = 0 (except at x = 0). This means that at t = 0 all load is carried
by the ¬‚uid and not by the solid. The aggregate modulus satis¬es K =
105 [Nm’2 ]; the permeability is given by H = 10’14 [m4 N’1 s’1 ].
(a) Adjust demo_fem1dcd to calculate the pressure as a function of
time for the above given problem. Choose v = 0 and take as a time
step t = 2500 [s]. Use θ = 1 for the time integration scheme.
(b) Determine the ¬‚uid pressure near the contact surface with the indenter
and plot this pressure as a function of time.
16 Solution of the three-dimensional
convection-diffusion equation by
means of the Finite Element Method

16.1 Introduction

The two- and three-dimensional convection-diffusion equation plays an important
role in many applications in biomedical engineering. One typical example from
recent research is the analysis of the effectiveness of different types of bioreactors
for tissue engineering. Tissue engineering is a rapidly evolving interdisciplinary
research area aiming at the replacement or restoration of diseased or damaged tis-
sue. In many cases devices made of arti¬cial materials are only capable of partially
restoring the original function of native tissues, and may not last for the full life-
time of a patient. In addition, there is no arti¬cial replacement for a large number
of tissues and organs. In tissue engineering new, autologous tissues are grown out-
side the human body by seeding cultured cells on scaffolds and further developed
in a bioreactor for later implantation. The tissue proliferation and differentiation
process is strongly affected by mechanical stimuli and transport of oxygen, min-
erals, nutrients and growth factors. To optimize bioreactor systems it is necessary
to analyse how these systems behave. The convection-diffusion equation plays an
important role in this kind of simulating analysis.
Fig. 16.1 shows two different bioreactor con¬gurations, which both have been
used in the past to tissue engineer articular cartilage. The work was especially
focussed on glucose, oxygen and lactate, because these metabolites play a major
role in the chondrocyte biosynthesis and survival. Questions ranged from: ˜Does
signi¬cant nutrient depletion occur at the high cell concentrations required for
chondrogenesis?™ to ˜Do increasing transport limitations due to matrix accumula-
tion signi¬cantly affect metabolite distributions?™ Fig. 16.2 shows a typical result
of the calculated oxygen distributions in the two bioreactor con¬gurations.
This chapter explains the discretization of the convection-diffusion equation
in two or three dimensions. First the diffusion equation is discussed, thereafter
the convection-diffusion equation is elaborated. The spatial discretization of the
weighting function is based on the Galerkin method.
The three-dimensional convection-diffusion equation

construct medium



Figure 16.1
Culture con¬gurations. (a) Petri dish (b) Compression set-up. Adapted from [16].

Oxygen [mMol]







Figure 16.2
Oxygen distribution in the static case at 48 h. Because of axisymmetry only the right half of the
domain cross section is shown. The construct position is indicated with a white line. (a) Petri dish
(b) Compression set-up. Adapted from [16].

16.2 Diffusion equation

Consider a two- or three-dimensional domain with boundary . As in the one-
dimensional model problem, the boundary is split into a part u along which the
essential boundary conditions are speci¬ed, and a part p along which the natural
boundary conditions may be speci¬ed. The generic form of the diffusion equation
is given by

∇ · ( c∇u) + f = 0, (16.1)
279 16.3 Divergence theorem and integration by parts

where c denotes the diffusion coef¬cient and f a source term. A more general form
of Eq. (16.1) is obtained by replacing the scalar c with a second-order tensor:

∇ · ( C · ∇u) + f = 0 . (16.2)

However, currently attention is restricted to Eq. (16.1).

The essential boundary conditions along read:

u = U at u, (16.3)

while the natural boundary conditions along are given by

n · c∇u = P at p, (16.4)

with n the unit outward normal vector to the boundary .

16.3 Divergence theorem and integration by parts

, and φ a
Let n be the unit outward normal to the boundary of the domain
suf¬ciently smooth function on , then

∇φ d = nφ d . (16.5)

If the function φ is replaced by a vector it can easily be derived:

∇ ·φd = n·φd . (16.6)

Eq. (16.6) is known as the divergence theorem. For a proof of these equations,
see for example Adams [1].
Let both φ and ψ be suf¬ciently smooth functions on , then

( ∇φ) ψ d = nφψ d ’ φ ∇ψ d . (16.7)

This is called integration by parts. To prove this we must integrate the product rule
of differentiation:

∇( φψ) = ( ∇φ) ψ + φ ∇ψ, (16.8)

to obtain

∇( φψ) d = ( ∇φ) ψ d + φ ∇ψ d . (16.9)
The three-dimensional convection-diffusion equation

Subsequently, use the divergence theorem to convert the left-hand side into the
boundary integral:

∇( φψ) d = nφψ d . (16.10)

This yields the desired result.

16.4 Weak form

Following the same steps as in Chapter 14, the differential equation Eq. (16.1) is
multiplied with a weighting function w and integrated over the domain :

w ∇ · ( c∇u) + f = 0,
d for all w. (16.11)

Next the integration by parts rule according to Eq. (16.7) is used:

w n · c∇u d ’ ∇w · ( c∇u) d + = 0.
wf d (16.12)

The boundary integral can be split into two parts, depending on the essential and
natural boundary conditions:

w n · c∇u d = w n · c∇u d + wP d , (16.13)
u p

where c∇u · n = P at p is used. It will be clear that, similar to the derivations
in Chapter 14, the ¬rst integral on the right-hand side of Eq. (16.13) is unknown,
while the second integral offers the possibility to incorporate the natural bound-
ary conditions. For the time being we keep both integrals together (to limit the
complexity of the equations and rewrite Eq. (16.12) according to:

∇w · ( c∇u) d = + w n · c∇u d .
wf d (16.14)

16.5 Galerkin discretization

Step 1 Introduce a mesh by splitting the domain into a number of non-
overlapping elements e . In a two-dimensional con¬guration the elements typ-
ically have either a triangular (in this case the mesh is sometimes referred to as a
triangulation) or a quadrilateral shape. A typical example of triangulation is given
in Fig. 16.3. Each triangle corresponds to an element.
281 16.5 Galerkin discretization

Figure 16.3
Example of a two-dimensional ¬nite element mesh using triangular elements. One element has
been highlighted.

The integration over the domain can be performed by a summation of the
integrals over each element:
Nel Nel
∇w · ( c∇u) d = + wn · c∇u d
wf d . (16.15)
e e e
e=1 e=1

The boundary part e denotes the intersection of element e with the boundary
, hence e = e © . Clearly, not each element will have an intersection with .

Step 2 If a Cartesian coordinate system is used for two-dimensional problems
(extension to the three-dimensional case is straightforward), the inner product
∇w · ( c∇u) yields
‚w ‚w ‚u ‚u
∇w · ( c∇u) = ex + ey · c ex + e y
‚x ‚y ‚x ‚y
‚w ‚u ‚w ‚u
=c + . (16.16)


. 56
( 67 .)