. 53
( 67 .)


(e) Show that the source term leads to a right-hand side column on
element level given by

he 1
fe = ’ .
∼ 2

(f) After discretization the resulting set of equations is expressed as

K u = f.
∼ ∼
261 Exercises

De¬ne the solution array u and derive the coef¬cient matrix K and the

array f for the two element mesh depicted in the ¬gure.

(g) Determine the solution array u.∼
14.6 Using an isoparametric formulation, the shape functions are de¬ned with
respect to a local coordinate system ’1 ¤ ξ ¤ 1. Within an element the
unknown uh is written as
uh = Ni ( ξ ) ui .

(a) What are the above shape functions with respect to the local coordi-
nate system if a linear or quadratic interpolation is used?
(b) How is the derivative of the shape functions
Assume a quadratic shape function, and let x1 = 0, x2 = 1, and
x3 = 3. Compute
for i = 1, 2, 3.
(d) Compute from array h given by

3 dN
h= ∼
x dx,


the ¬rst component using the same quadratic shape functions as
14.7 Consider in the code mlfem_nac the directory oneD. The one-
dimensional ¬nite element program fem1d solves the diffusion problem:
+ f = 0,
dx dx
on a domain a ¤ x ¤ b subject to given boundary conditions for a cer-
tain problem. The input data for this program are speci¬ed in the m-¬le
demo_fem1d, along with the post-processing statements.
(a) Modify the m-¬le demo_fem1d such that ¬ve linear elements are
used to solve the above differential equation, using the boundary

u = 0 at x = 0,
Numerical solution of one-dimensional diffusion equation

p=c = 1 at x = 1,

with c = 1 and f = 0. Compute the FEM solution (array sol) and
compare the nodal values with the exact solution.
(b) Use the array pos to extract the nodal solutions within the second
element from the array sol.
(c) Determine the solution for the third node using the array dest.
14.8 One of the major problems in the tissue engineering of cartilage is to make
thick constructs. Nutrients coming from the surrounding medium have to
reach the cells in the middle of the construct by diffusion, but cells close
to the edge of the construct may consume so much nutrients that there is
nothing left for cells in the middle. Consider the experimental set-up as
given in the ¬gure below, representing a schematic drawing of a bioreactor
to culture articular cartilage.
The construct is ¬xed between two highly permeable membranes allowing
free contact with the culture medium. The thickness t is a trade-off between
the diffusion coef¬cient c, the consumption rate of the cells f (both cannot
be in¬‚uenced) and the amount of the molecules that can be supplied via the
medium, which is usually bound to a maximum.
The diffusion problem can be considered as a one-dimensional problem.
The current analysis is meant to determine the glucose concentration
u( x), which is an essential nutrient for the cells. The following proper-
ties are given: c = 9.2 —10’6 [cm2 s’1 ] and f = 56 — 10’7 [Mol hour’1
cm’3 ]. Low glucose medium is used, concentration = 5 — 10’3 [Mol
litre’1 ].
(a) Assuming the consumption rate f is constant and the medium is
refreshed continuously, a stable glucose concentration as a func-
tion of the location in the construct will be reached after a while.
Give the differential equation and boundary conditions describing this
(b) Give the analytical solution of this problem by integrating the
differential equation.
(c) Adjust demo_fem1d to solve the problem with the ¬nite element
method. Solve the problem for t = 1, 2, 3, 4 and 5 [mm].
(d) At which thickness do the cells appear to die in the middle of the
14.9 To determine the material properties of a skeletal muscle a uniaxial tensile
test is performed as shown in the ¬gure. The muscle is clamped on one side
and a force F = 10 [N] is applied on the other side. The muscle has a total
263 Exercises

Medium Medium


glucose concentration
u(x )


length of = 12 [cm]. The muscle has a circular cross section. The radius
of the cross section can be approximated by:
r = a1 sin3 ( a2 x + a3 ) ,


with a1 = 1.6 [cm], a2 = 0.15 [cm’1 ] and a3 = 0.8 [-]. The estimated
Young™s modulus is E = 105 [N m’2 ].
(a) Give the differential equation for the axial displacement u( x) and
boundary conditions that describe the current problem.
(b) Determine the displacement ¬eld by adjusting the ¬le demo_fem1d.
15 Solution of the one-dimensional
convection-diffusion equation by
means of the Finite Element Method
15.1 Introduction

This chapter extends the formulation of the previous chapter for the one-
dimensional diffusion equation to the time-dependent convection-diffusion equa-
tion. Although a good functioning of the human body relies on maintaining a
homeostasis or equilibrium in the physiological state of the tissues and organs,
it is a dynamic equilibrium. This means that all processes have to respond to
changing inputs, which are caused by changes of the environment. The diffusion
processes taking place in the body are not constant, but instationary, so time has
to be included as an independent variable in the diffusion equation. Thus, the
instationary diffusion equation becomes a partial differential equation.
Convection is the process whereby heat or particles are transported by air or
¬‚uid moving from one point to another point. Diffusion could be seen as a process
of transport through immobilized ¬‚uid or air. When the ¬‚uid itself moves, particles
in that ¬‚uid are dragged along. This is called convection and also plays a major
role in biomechanics. An example is the loss of heat because moving air is passing
the body. The air next to the body is heated by conduction, moves away and carries
off the heat just taken from the body. Another example is a drug that is released
at some spot in the circulation and is transported away from that spot by means of
the blood ¬‚ow. In larger blood vessels the prime mechanism of transportation is

15.2 The convection-diffusion equation

Assuming that the source term f = 0, the unsteady one-dimensional convection-
diffusion equation can be written as
‚u ‚u ‚ ‚u
+v = c , (15.1)
‚t ‚x ‚x ‚x
with u a function of both position x and time t:
u = u(x, t) . (15.2)
265 15.2 The convection-diffusion equation

The convective velocity is denoted by v and c is the diffusion coef¬cient. Both v
and c are assumed to be constant in the present chapter. Compared to the dif-
fusion equation of the previous chapter, two terms have been added: the time
dependency term ‚u/‚t (inertia term) and the convective term v‚u/‚x. The
=[ a, b],
convection-diffusion equation holds within a given spatial domain
i.e. with boundaries at x = a and x = b, as well as within a time domain, say
S = [ 0, T], and is assumed to be subject to the boundary conditions:

(located at x = a) ,
u = U at (15.3)

= P at p (located at x = b) .
c (15.4)
Furthermore, one initial condition on u must be speci¬ed, say at t = 0:

u( x, t = 0) = uini ( x) in . (15.5)

Under certain conditions the transient character of the solution of the unsteady
convection-diffusion equation vanishes. In that case we deal with a steady
convection-diffusion problem, described by the differential equation
d du
= c . (15.6)
dx dx dx

Example 15.1 A typical example of a solution of the steady convection-diffusion problem is
represented in Fig. 15.1. In this example the domain spans 0 ¤ x ¤ 1, while at
x = 0 the solution is set to 0 and at x = 1 the value u = 1 is imposed. That means
that in this example two essential boundary conditions are used and that there is
no natural boundary condition prescribed. The analytical solution for v = 0 in this
case is
1 v
v (1 ’ e ) .
u= cx
1 ’ ec
Clearly, without any convection, i.e. v = 0, a spatially linear distribution of u
results, while with increasing convective velocity v a boundary layer develops at
x = 1.
The convection-diffusion equation may be written in a dimensionless form by
introducing an appropriate length scale L, a time scale and a reference solution
U, for example the reference temperature or concentration. Then, the dimension-
less solution u— , the dimensionless coordinate x— and the dimensionless time t—
are de¬ned as
u x t
u— = x— = t— = . (15.7)
The one-dimensional convection-diffusion equation













. 53
( 67 .)