· are easily computed. However, the shape function derivatives with respect

to the global coordinates x and y are needed. For this purpose the concept of

isoparametric elements is used.

For isoparametric elements the global coordinates within an element are

interpolated based on the nodal coordinates using the shape functions

= N T ( ξ , ·) xe , = N T ( ξ , ·) ye ,

y| (16.44)

x| ∼

∼ ∼

e e ∼

where xe and ye contain the nodal x- and y-coordinates, respectively. These equa-

∼ ∼

tions re¬‚ect the transformation from the local coordinates ( ξ , ·) to the global

coordinates ( x, y). The derivatives of Ni with respect to the Cartesian coordinates

x and y can be evaluated with the aid of the chain rule:

The three-dimensional convection-diffusion equation

286

1

0.8

0.6

N1

1

0.4

0.2

0 ·

0

“1

“0.5

0 “1

0.5

ξ 1

Figure 16.5

Shape function associated with node 1.

‚Ni ‚Ni ‚ξ ‚Ni ‚·

= +

‚ξ ‚x ‚· ‚x

dx

‚Ni ‚Ni ‚ξ ‚Ni ‚·

= + . (16.45)

‚ξ ‚y ‚· ‚y

dy

These relations can be rewritten in the following matrix form:

⎡ ¤⎡ ¤⎡ ¤

‚Ni ‚ξ ‚· ‚Ni

⎢ ‚x ⎥ ⎢ ‚x ‚x ⎥ ⎢ ‚ξ ⎥

⎢ ⎥⎢ ⎥⎢ ⎥.

⎣ ‚Ni ¦ = ⎣ ‚ξ (16.46)

‚· ¦ ⎣ ‚Ni ¦

‚y ‚y ‚y ‚·

The derivatives ‚Ni /‚ξ and ‚Ni /‚· are readily available, but the terms in the

matrix cannot be directly computed since the explicit expressions ξ = ξ ( x, y) and

· = ·( x, y) are not known. However, due to the isoparametric formulation the

inverse relations are known, so the following matrix can be calculated easily:

⎡ ¤

‚x ‚y

⎢ ‚ξ ‚ξ ⎥

x,ξ = ⎢ ⎥, (16.47)

⎣ ‚x ‚y ¦

‚· ‚·

287 16.7 Isoparametric elements and numerical integration

with components:

‚x ‚N T

= ∼

x

‚ξ ∼e

‚ξ

‚x ‚N T

= ∼

x

‚· ∼e

‚·

(16.48)

‚y ‚N T

= ∼

y

‚ξ ‚ξ ∼e

‚y ‚N T

= ∼

y.

‚· ‚· ∼e

Matrix (16.47) is the inverse of the matrix in Eq. (16.46) (this can be checked by

multiplying the two matrices, which gives the unit matrix). Accordingly

⎡ ¤ ⎡ ¤

‚·

‚ξ ‚y ‚y

’ ‚ξ

⎣ ‚x ‚x ¦ = ( x,ξ )’1 = 1 ⎣ ‚· ¦ (16.49)

‚·

‚ξ ‚x ‚x

’

j

‚· ‚ξ

‚y ‚y

where

‚x ‚y ‚x ‚y

j = det( x,ξ ) = ’ . (16.50)

‚ξ ‚· ‚· ‚ξ

It is an elaborate process and usually not possible to analytically compute the inte-

grals in the expressions for the element matrices or arrays, so generally they are

approximated by numerical integration. Each of the components of the matrices,

such as K e etc., that need to be computed consists of integrals of a given function,

say g( x, y), over the domain of the element e . These may be transformed to an

integral over the unit square ’1 ¤ ξ ¤ 1, ’1 ¤ · ¤ 1, according to

1 1

= f ( x( ξ , ·) , y( ξ , ·) ) j( ξ , ·) dξ d·,

f ( x, y) d (16.51)

’1 ’1

e

with j( ξ , ·) de¬ned by Eq. (16.50).

The integral over the unit square may be approximated by a numerical

integration (quadrature) rule, giving

nint

1 1

f ( ξ , ·) j( ξ , ·) dξ d· ≈ f ( ξi , ·i ) j( ξi , ·i ) W( ξi , ·i ) . (16.52)

’1 ’1 i=1

For example, in case of a two-by-two Gaussian integration rule the location of the

integration points have ξ , ·-coordinates and associated weights:

’1 ’1

ξ1 = √ , ·1 = √ , W1 = 1

3 3

’1

1

ξ2 = √ , ·2 = √ , W2 = 1

3 3

(16.53)

The three-dimensional convection-diffusion equation

288

1 1

ξ3 = √ , ·3 = √ , W3 = 1

3 3

’1 1

ξ4 = √ , ·4 = √ , W4 = 1.

3 3

16.8 Example

One of the treatments for coronary occlusions that may lead to an infarct is to put a

stent at the location of the occlusion. One of the problems with this intervention is

that the blood vessels often occlude again quite soon after the stent is placed. One

of the solutions may be to design a stent that gradually releases drugs to prevent

such an occlusion from occurring again, see Fig. 16.6. How these drugs propagate

through the vascular tree is a convection-diffusion problem.

Consider the domain as sketched in Fig. 16.7. It represents a section of a long

channel. For reasons of simplicity the three-dimensional problem is modelled as a

two-dimensional problem: the relevant ¬elds in the con¬guration are assumed to

be independent of the coordinate perpendicular to the xy-plane. A so-called New-

tonian ¬‚uid (modelling blood in a ¬rst approximation) ¬‚ows through the channel

as indicated in the ¬gure.

Let u denote the concentration of a certain drug. Along the entrance of the

domain the drug concentration is zero. Along the small part of the wall indicated

Figure 16.6

Schematic of a stent in a blood vessel.

y

2

1

0

“1

“2

“2 0 2 4 6 8 10 12 x

Figure 16.7

Speci¬cation of the convection-diffusion problem.

289 16.8 Example

2

1.5

C2

1

C1

0.5

y

0

’0.5

’1

’1 0 1 2 3 4 5 6 7 8 9

x

Figure 16.8

Computational domain of the convection-diffusion problem.

with a thick line, the drug concentration is prescribed, say u = 1. The drug diffuses

into the liquid with a diffusion constant c, but is also convected by the ¬‚uid. The

aim is to compute the concentration pro¬le in the two-dimensional channel for a

number of ¬‚uid velocities. The computational domain is indicated by the dashed

line, and is further outlined in Fig. 16.8. Because of symmetry only the top half of

the vessel is modelled.

For stationary ¬‚ow conditions, the velocity ¬eld v is described by means of a

parabolic pro¬le (Poisseuille ¬‚ow) according to

v = a(1 ’ y2 ) ex .

As mentioned before, along boundary C1 the ¬‚uid ¬‚ows into the domain with

a concentration u = 0, while along boundary C2 the concentration u = 1

is prescribed. Along the remaining parts of the boundary the natural boundary

condition:

n · c∇u = 0,

is imposed. This means that the top wall is impenetrable for the drug, while this

condition must also be enforced along the symmetry line y = 0. Speci¬cation

of this condition on the out¬‚ow boundary is somewhat disputable, but dif¬cult to

avoid, because only a small part of the circulation system is modelled. By choos-

ing the out¬‚ow boundary far away from the source of the drug the in¬‚uence of

this boundary condition is small.

The corresponding mesh is shown in Fig. 16.9. The problem is discretized using

bi-quadratic elements.

The steady convection diffusion problem is solved, for c = 1 and a sequence

0,1,10,25,100 of parameter a. Clearly, with increasing a the velocity in the x-

direction increases proportionally, hence convection becomes increasingly impor-

tant. For increasing a contours of constant u are depicted in Fig. 16.10. In all

The three-dimensional convection-diffusion equation

290

2

1.5

C2

1

C1

0.5

y

0

’0.5

’1

’1 0 1 2 3 4 5 6 7 8 9

x

Figure 16.9

Mesh of the convection-diffusion problem.