element level given by

he 1

fe = ’ .

1

∼ 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

n

uh = Ni ( ξ ) ui .

i=1

(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

dNi

dx

obtained?

Assume a quadratic shape function, and let x1 = 0, x2 = 1, and

(c)

x3 = 3. Compute

dNi

for i = 1, 2, 3.

dx

(d) Compute from array h given by

∼

3 dN

h= ∼

x dx,

∼

dx

0

the ¬rst component using the same quadratic shape functions as

above.

14.7 Consider in the code mlfem_nac the directory oneD. The one-

dimensional ¬nite element program fem1d solves the diffusion problem:

du

d

+ f = 0,

c

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

conditions

u = 0 at x = 0,

and

Numerical solution of one-dimensional diffusion equation

262

du

p=c = 1 at x = 1,

dx

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

process.

(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

construct?

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

construct

Medium Medium

t

glucose concentration

u(x )

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 ) ,

F

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

convection.

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)

u

and

‚u

= P at p (located at x = b) .

c (15.4)

‚x

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

du

= c . (15.6)

v

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)

U L

The one-dimensional convection-diffusion equation

266

1

0.9

0.8

0.7

0.6

v

0.5

u

0.4

0.3

0.2

0.1

0