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

filter

supporting

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

276

that the pressure in the ¬‚uid is described by the instationary diffusion

equation:

‚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

278

construct medium

(a)

medium

construct

(b)

Figure 16.1

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

Oxygen [mMol]

0.18

0.14

(a)

0.1

0.06

0.02

(b)

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 = U at u, (16.3)

while the natural boundary conditions along are given by

p

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

280

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)