a=0

(a)

a=1

(b)

a = 10

(c)

a = 25

(d)

a = 100

(e)

Figure 16.10

Contour lines of constant u values for a range of values of a. With increasing a the effect of

convection increases.

291 Exercises

cases 10 contour lines are shown ranging from 0.1 to 1 with increments of 0.1.

The effect of an increasing velocity is clearly demonstrated.

Exercises

16.1 The weak form of the two-dimensional convection-diffusion equation is

given by

‚u

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

w d

‚t

= + wn · c∇u = d .

wf d

After discretization, the element inertia matrix is de¬ned as

Me = N NT d ,

∼∼

e

while the element matrix related to the convective part is given by

dN T dN T

Ce = + vy ∼

∼

N vx d.

∼

dx dy

e

Consider the element as depicted in the ¬gure below.

y

4 3

x

1 2

This is a bi-linear element. The element spans the spatial domain ’2 ¤

x ¤ 2 and 2 ¤ y ¤ 2.

(a) Compute the element inertia matrix M e using a 2—2 Gauss integration

rule. It is recommended to use MATLAB for this computation.

(b) Suppose that the location of the integration points coincides with the

nodes of the element. What is the element inertia matrix M e in this

case?

(c) Compute the matrix Ce if vx = 1 and vy = 0, using a 2 — 2 Gauss

integration rule.

(d) Suppose that along the edge of the element located at x = 2 a constant

¬‚ux q = n · c∇u is prescribed. Then compute the column

The three-dimensional convection-diffusion equation

292

Nq d ,

∼

e

for this edge, contributing to the element column f e .

∼

16.2 Consider the mesh depicted in the ¬gure below. The solution vector u con-

∼

tains the nodal solutions in the sequence de¬ned by the node numbers,

hence:

uT = [ u1 u2 · · · u16 ] .

∼

10 12 14 16

3 6 9

5 7 9 15

2 5 8

13

2 4 8

1 4 7

6 11

1 3

Let the element topology array of the third element be given by

top(3,:)=[5 7 12 10 1 1].

(a) Is this topology array unique?

(b) What is the pos array for this element?

(c) De¬ne the dest array.

(d) The solution u is stored in the array sol. How are the nodal solutions

∼

of element 8 extracted from sol using the array pos?

(e) Suppose that the array nodes contains the node numbers of the

left edge, say nodes=[1 2 5 10]. How are the nodal solutions

of these nodes extracted from the solution array sol via the array

dest?

(f) Consider the third element whose element topology array is given by

top(3,:)=[5 7 12 10 1 1]. If the nodal solution array ue for ∼

this element is given by:

uT = [ 1 4 3 7] ,

∼e

compute the solution u for ξ = 1 and · = 3 .

4 4

(i) Suppose the element topology array is given by

top(3,:)=[12 10 5 7 1 1]. What is the corresponding ele-

ment nodal solution array ue ?

∼

293 Exercises

16.3 Consider the square domain 0 ¤ x ¤ 1, 0 ¤ y ¤ 1. Along the line x = 0

the boundary condition u = 0 is imposed, and along the line x = 1 the

boundary condition u = 1 is prescribed. Along the other boundaries the

natural boundary condition n · c∇u = 0 is speci¬ed. Consider the steady

convection diffusion problem on this domain.

(a) Suppose that the convection-diffusion problem represents the tem-

perature equation. What physical meaning does the above natural

boundary condition have?

(b) Adjust the m-¬le demo_cd such that the above problem is solved

using 5 — 5 linear elements. Select the diffusion constant c = 1 and

the convective velocity v = ex (Remark: see the element m-¬le elcd

and the associated m-¬le elcd_a). Specify the structured array mat

according to this.

(c) Extract the solution along the line y = 0. Hint, ¬rst select the nodes

using usercurves containing the nodes along all the curves that

have been de¬ned for the mesh generator crmesh. Then use dest

to extract the relevant solution components.

(d) Compare the solution along y = 0 with the one-dimensional solution,

using the program fem1dcd.

16.4 Consider the domain as sketched in the ¬gure. Along the boundary 1 the

speci¬cation u = 1 is chosen, while along 7 u = 0 is chosen as a boundary

condition. Along all other boundaries n · c∇u = 0 is imposed. Choose the

convective velocity v such that:

y<0: v= 0

y≥0: v = y ex

Furthermore, choose c = 1 throughout the domain. Model this problem

using femlin_cd. Approximately beyond which value of x along the

boundary 6 holds u > 0.8?

3

“6

2

“7 “5

1

y

“3

0

“4 “2

“1

’1

’2

0 2 4 6 8 10

x

The three-dimensional convection-diffusion equation

294

16.5 If the ¬‚ux n · c∇u = P along is prescribed, the integral (element level):

p

NP d ,

∼

e

needs to be computed. It may be transformed to the local coordinate system

ξ ∈[ ’1 1], such that:

‚x

1

NP d = NP dξ ,

‚ξ

∼ ∼

’1

e

where it is assumed that e is oriented in x-direction.

If linear elements are chosen, the shape functions along the boundary are

given by:

2 (1 ’ ξ )

1

N= .

2 (1 + ξ )

1

∼

Use an isoparametric element and let L denote the length of e, then

demonstrate that for constant P:

LP 1

NP d = .

∼

1

2

e

17 Shape functions and numerical

integration

17.1 Introduction

In the previous chapter the shape functions Ni have hardly been discussed in any

detail. The key purpose of this chapter is ¬rst to introduce isoparametric shape

functions, and second to outline numerical integration of the integrals appearing

in the element coef¬cient matrices and element column. Before this can be done

it is useful to understand the minimum requirements to be imposed on the shape