c=1

u( x = 0) = 0

u( x = 1) = 1

The convective velocity v will be varied. For v = 0, the diffusion limit, the solution

is obvious: u varies linearly in x from u = 0 at x = 0 to u = 1 at x = 1.

Figs. 15.3(a) to (d) show the solution for v = 1, 10, 25 and 100, respectively,

using a uniform element distribution with ten linear elements. For v = 1 and v =

10 the approximate solution uh (solid line) closely (but not exactly) follows the

exact solution (dashed line). However, for v = 25 the numerical solution starts to

demonstrate an oscillatory behaviour that is more prominent for v = 100. Careful

analysis of the discrete set of equations shows that the so-called element Peclet

number governs this oscillatory behaviour. The element Peclet number is de¬ned

as

vh

Peh = ,

2c

271 15.4 Spatial discretization

1 1

0.9 0.9

0.8 0.8

0.7 0.7

0.6 0.6

u 0.5 u 0.5

0.4 0.4

0.3 0.3

0.2 0.2

0.1 0.1

0 0

0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

x x

(a) (b)

1

1.2

0.8

1

0.6

0.8 0.4

0.6 0.2

u

u

0

0.4

“0.2

0.2

“0.4

0 “0.6

“0.2 “0.8

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

x

x

(c) (d)

Figure 15.3

Solution of the steady convection-diffusion equation for v = 1, 10, 25 and v = 100, respectively

using ten linear elements; solid line: approximate solution uh , dashed line: exact solution u.

where h is the element length. Above a certain critical value of Peh the solu-

tion behaves in an oscillatory fashion. To reduce possible oscillations the element

Peclet number should be reduced. For ¬xed v and c this can only be achieved by

reducing the element size h. For example, doubling the number of elements from

10 to 20 eliminates the oscillations at v = 25, see Fig. 15.4.

The oscillations that appear in the numerical solution of the steady convection-

diffusion equation may be examined as follows. Consider a domain that is

subdivided in two linear elements, each having a length equal to h. At one end

of the domain the solution is ¬xed to u = 0, while at the other end the solution

is set to u = 1, or any other arbitrary non-zero value. For constant v and c the

governing differential equation may be rewritten as

v du d2 u

’ 2 = 0.

c dx dx

The one-dimensional convection-diffusion equation

272

1

0.9

0.8

0.7

0.6

0.5

u

0.4

0.3

0.2

0.1

0

0 0.2 0.4 0.6 0.8 1

x

Figure 15.4

Solution of the steady convection-diffusion problem using 20 elements at v = 25.

The set of equations that results after discretization is, as usual:

Ku=f .

∼ ∼

If only two linear elements of equal length h are used the coef¬cient matrix K may

be written as

⎡ ¤ ⎡ ¤

’1 1 0 1 ’1 0

v⎢ ⎥ 1⎢ ⎥

K= ⎣ ’1 0 1 ¦ + ⎣ ’1 2 ’1 ¦ ,

2c h

0 ’1 1 0 ’1 1

where the ¬rst, asymmetric, part corresponds to the convective term and the sec-

ond, symmetric, part to the diffusion term. In the absence of a source term the

second component of f is zero. Let u1 and u3 be located at the ends of the domain

∼

such that u1 = 0 and u3 = 1, then u2 is obtained from

vh

2u2 = 1 ’

.

2c

An oscillation becomes manifest if u2 < 0. To avoid this, the element Peclet

number should be smaller than one:

vh

Peh = < 1.

2c

Consequently, at a given convective velocity v and diffusion constant c, the mesh

size h can be chosen such that an oscillation-free solution results. In particular

for large values of v/c this may result in very ¬ne meshes. To avoid the use of

273 Exercises

t = 0.01 : 0.05 : 0.26

1

0.9

0.8

0.7

0.6

0.5

u

0.4

0.3

0.2

0.1

0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

x

Figure 15.5

Solution of the unsteady convection-diffusion problem using 20 elements at v = 10.

very ¬ne meshes, an alternative, stabilized formulation has been developed: the

so-called SUPG (Streamline-Upwind/Petrov“Galerkin) formulation. A discussion

of this method however, is beyond the scope of the present book.

Example 15.3 Let us consider the instationary convection-diffusion problem. For this prob-

lem the same outline as in the previous example for v = 10 is chosen and

we use a uniform distribution of 20 linear elements. The initial condition is

u( x, t = 0) = 0 throughout the domain. At the ¬rst time step the boundary

condition u( x = 1, t) = 1 is imposed. The unsteady solution is obtained using a

time step of t = 0.01, while θ = 0.5 is selected for the θ-scheme. Fig. 15.5

shows the time-evolution of the solution towards the steady state value (denoted

by the dashed line) for v = 10.

Exercises

15.1 Consider the domain = [ 0 1]. On the domain the one-dimensional steady

convection-diffusion equation:

d2 u

du

=c 2

v

dx dx

The one-dimensional convection-diffusion equation

274

holds. As boundary conditions, at x = 0, u = 0 and at x = 1, u = 1 are

speci¬ed.

(a) Prove that the exact solution is given by

1 v

u= ( 1 ’ e c x) .

v

1’e c

(b)Verify this by means of the script demo_fem1dcd, to solve the one-

dimensional convection-diffusion problem, which can be found in the

directory oned of the programme library mlfem_nac. Use ¬ve ele-

ments and select c = 1, while v is varied. Choose v = 0, v = 1,

v = 10 and v = 20. Explain the results.

(c) According to Section 15.4, the solution is expected to be oscillation

free if the element Peclet number is smaller than 1:

ah

Peh = < 1.

2c

Verify that this is indeed the case.

15.2 Investigate the unsteady convection-diffusion problem:

‚u ‚u ‚ ‚u

+v = c ,

‚x ‚x ‚x

dt

= [ 0 1] subject to the initial condition:

on the domain

uini ( x, t = 0) = 0,

inside the domain and the boundary conditions:

u = 0 at x = 0, u = 1 at x = 1.

The θ-scheme for time integration is applied. Modify the m-¬le

demo_fem1dcd accordingly. Use ten linear elements.

(a) Choose v = c = 1 and solve the problem with different values of

θ. Use θ = 0.5, θ = 0.4 and θ = 0.25. For each problem start

with a time step of 0.01 and increase the time step with 0.01 until a

maximum of 0.05. Describe what happens with the solution.

(b) In the steady state case, what is the maximum value of the convective

velocity v such that the solution is oscillation free for c = 1?

(c) Does the numerical solution remain oscillation free in the unsteady

case for θ = 0.5 and t = 0.001, 0.01, 0.1? What happens?

(d) What happens at t = 0.001, if the convective velocity is reduced?

15.3 Investigate the unsteady convection problem:

‚u ‚ ‚u

‚u

+v = c

‚t ‚x ‚x ‚x

275 Exercises

on the spatial domain = [ 0 1] and a temporal domain that spans t =

[ 0 0.5]. At x = 0 the boundary condition u = 0 is prescribed, while at x =

1, c du/dx = 0 is selected. The convection dominated case is investigated,

with v = 1 and c = 0.01. The problem is solved using the θ-scheme with

θ = 0.5, a time step t = 0.05 and 40 linear elements.

(a) First, the initial condition:

u(x, t = 0) = sin(2π x)

is considered. Adapt the demo_fem1dcd and the FEM program

accordingly. In particular, make sure that the initial condition is han-

dled properly in fem1dcd. (Hint: the initial condition is speci¬ed in

sol(:,1).) Solve this problem. Plot the solution for all time steps.

(b) Second, let the initial condition be given by

u( x, t = 0) = 0 for x < 0.25 and x > 0.5,