<<

. 58
( 67 .)



>>

The shape function derivatives with respect to the local coordinates ξ and
· 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.

<<

. 58
( 67 .)



>>