Step 3 Introduce a discretization for both the weighting function w and the

unknown u, so within each element

n

wh | = Ni ( x, y) we,i = N T ( x, y) we (16.17)

∼ ∼

e

i=1

n

uh | = Ni ( x, y) ue,i = N T ( x, y) ue . (16.18)

∼

∼

e

i=1

The three-dimensional convection-diffusion equation

282

Step 4 Substitution of this discretization into Eq. (16.16) yields

‚wh ‚uh ‚wh ‚uh

∇wh · c∇uh = c +

‚x ‚x ‚y ‚y

‚N T we ‚N T ue ‚N T w ‚N T u

+ ∼ ∼ e ∼ ∼e

=c ∼∼

∼∼

‚x ‚x ‚y ‚y

T ‚N ‚N ‚N ‚N T

T

= c we +∼∼

∼ ∼

u. (16.19)

‚y ‚y ∼e

‚x ‚x

∼

Step 5 Using the discretization in Eq. (16.15) yields

Nel

‚N ‚N T ‚N ‚N T

+∼∼

wT ∼ ∼

c d ue

∼e

‚x ‚x ‚y ‚y ∼

e

e=1

Nel

= + wT N n · c∇u d

wT Nf d . (16.20)

∼e ∼e

∼ ∼

e e

e=1

The element matrix is given by

‚N ‚N T ‚N ‚N T

Ke = +∼∼

∼ ∼

c d, (16.21)

‚x ‚x ‚y ‚y

e

and the element column by

f= + N n · c∇u d .

Nf d (16.22)

∼ ∼

∼e

e e

Using this notation, Eq. (16.20) may be written as

Nel Nel

ue =

wT K e wT f e . (16.23)

∼e ∼e

∼ ∼

e=1 e=1

Following a similar procedure as outlined in Chapter 14, this may be rearranged

into

w T K u = wT f , (16.24)

∼

∼ ∼ ∼

which should hold for all w. This ¬nally leads to:

∼

Ku = f . (16.25)

∼ ∼

Similar to the situation described in Section 14.6 the column u contains an

∼

unknown and a known part depending on the essential boundary conditions. In the

nodes, where essential boundary conditions are prescribed, the associated external

loads in the right-hand side integrals are unknown. However, the set equations can

be partitioned as has been done in Chapter 14 to arrive at a set that can be solved.

283 16.6 Convection-diffusion equation

16.6 Convection-diffusion equation

Assuming isotropic diffusion, the convection-diffusion equation is given by

‚u

+ v · ∇u = ∇ · ( c∇u) + f , (16.26)

‚t

with v the convective velocity. This equation should hold on the spatial domain

during a certain period of time, say S = [ 0, T]. Initial boundary conditions must

be speci¬ed:

u( x, t = 0) = uini ( x) in , (16.27)

as well as essential and natural boundary conditions:

u = U at (16.28)

u

n · c∇u = P at p. (16.29)

The weak form is obtained analogously to the procedure of Section 16.4, giving

‚u

+ wv · ∇u + ∇w · ( c∇u)

w d

‚t

= + wn · c∇u d .

wf d (16.30)

Spatial discretization is performed in a two-dimensional con¬guration by intro-

ducing

wh | = N T ( x, y) we , (16.31)

∼ ∼

e

and

uh | = N T ( x, y) ue . (16.32)

∼

∼

e

For a particular element e, the individual integrals of Eq. (16.30) can be

converted to:

‚u due

= wT M e ∼

w d , (16.33)

∼e

‚t dt

e

with:

Me = NNT d . (16.34)

∼∼

e

Further, by using

v = v x e x + vy e y , (16.35)

it can be written:

wv · ∇u d = wT C e u e , (16.36)

∼e ∼

e

The three-dimensional convection-diffusion equation

284

with

‚N T ‚N T

Ce = + vy ∼

∼

N vx d. (16.37)

‚x ‚y

∼

e

The remaining integrals from Eq. (16.30) follow from the previous section.

Therefore, Eq. (16.30) may be formulated as

Nel Nel

du

M e ∼e + ( Ce + K e ) ue =

wT wT f e , (16.38)

∼e ∼e

∼ ∼

dt

e=1 e=1

which must be satis¬ed for all admissible weighting values, hence after the

assembly process the following set of equations results:

du

+ ( C + K) u = f .

∼

M (16.39)

∼ ∼

dt

Application of the θ-scheme for temporal discretization in the time increment

[ tn , tn+1 ] results in:

1

+ θC + θK un+1

M ∼

t

1

un + f θ ,

=M ’( 1 ’ θ) C’( 1 ’ θ) K (16.40)

∼ ∼

t

t = tn+1 ’ tn and

with

f θ = θf n+1 + ( 1 ’ θ) f n . (16.41)

∼ ∼ ∼

16.7 Isoparametric elements and numerical integration

To compute the element matrices M e , Ce and K e and the element array f e the shape

∼

functions N and the shape function derivatives with respect to the coordinates

∼

x and y, need to be available. To de¬ne the shape functions it is convenient to

introduce a local coordinate system.

In particular for an arbitrarily shaped quadrilateral element, for instance as

depicted in Fig. 16.4(a) it is dif¬cult or impossible to de¬ne explicitly the shape

functions with respect to the global coordinates x and y. However, de¬ning the

shape functions in a unit square domain, as depicted in Fig. 16.4(b) is straight-

forward. The local coordinate system de¬ned in Fig. 16.4(b) is chosen such that

’1 ¤ ξ ¤ 1 and ’1 ¤ · ¤ 1. Node one, for instance, has local coordinates

ξ = ’1 and · = ’1. Remember that the shape function of a certain node has

to be one at the spatial location of this node and zero at the spatial location of all

other nodes. Using the local coordinate system allows an elegant de¬nition of the

shape functions of a four-node quadrilateral element:

285 16.7 Isoparametric elements and numerical integration

·

3 (“1, 1) (1, 1)

4 4 3

ξ

ey 1

2

2

1

ex (1, “1)

(“1, “1)

(a) Quadrilateral element with (b) Quadrilateral element with

respect to a global coordinate system respect to a local coordinate system

Figure 16.4

Quadrilateral element with respect to global and local coordinate systems.

1

N1 = (1 ’ ξ ) ( 1 ’ ·)

4

1

N2 = (1 + ξ ) ( 1 ’ ·)

4 (16.42)

1

N3 = (1 + ξ ) ( 1 + ·)

4

1

N4 = (1 ’ ξ ) ( 1 + ·) .

4

An element having these shape functions is called a bi-linear element. Along the

edges of the element the shape functions are linear with respect to either ξ or ·.

Within the element, however, the shape functions are bi-linear with respect to ξ

and ·. For instance:

1

N1 =

(1 ’ ξ ’ · + ξ ·) . (16.43)

4

Fig. 16.5 shows N1 visualized as a contour plot.