x

Figure 15.1

Solution to steady 1-D convection-diffusion problem.

Assuming v and c constant, the dimensionless form may be written as

1 ‚u— ‚u— ‚u—

‚

+ Pe — = — , (15.8)

Fo ‚t— ‚x—

‚x ‚x

with the Fourier number given by

c

Fo = , (15.9)

L2

and the Peclet number given by

vL

Pe =. (15.10)

c

The Peclet number re¬‚ects the relative importance of convection compared to dif-

fusion. It will be demonstrated that with increasing Peclet number the numerical

solution of the convection-diffusion problem becomes more dif¬cult.

In the next section we will start by studying the time discretization of an

instationary equation. After that in Section 15.4 the spatial discretization of the

convection-diffusion equation will be discussed.

15.3 Temporal discretization

Many algorithms have been developed for the temporal discretization of the

convection-diffusion equation. Only one method is discussed here: the so-called

267 15.3 Temporal discretization

θ-scheme. Getting ahead of the spatial discretization in the next section, which

will lead to a set of linear differential equations, the time discretization is best

illustrated with a set of ¬rst-order linear differential equations:

du

+ A u = f (t),

∼

(15.11)

∼ ∼

dt

where A is a constant matrix and f ( t) a given column. To solve this time-dependent

∼

problem the time domain is split into a ¬nite number of time increments. Then,

attention is focussed on the increment with tn < t < tn+1 = tn + t. Assuming

that the solution un at time tn is known, the unknown un+1 at time tn+1 has to be

determined. The θ-scheme approximates these differential equations by

un+1 ’ un

+ θA un+1 + (1 ’ θ) A un = θf n+1 + (1 ’ θ) f n .

∼ ∼

(15.12)

∼ ∼ ∼ ∼

t

For θ = 0 this scheme reduces to the Euler explicit or forward Euler scheme,

while for θ = 1 the Euler implicit or backward Euler scheme results. Both of

these schemes are ¬rst-order accurate, i.e. O( t). This means that the accuracy

of the solution is linearly related to the size of the time step t. The accuracy

improves when the time step becomes smaller. For θ = 0.5 the Crank“Nicholson

scheme results which is second-order accurate, i.e. O( t2 ).

To illustrate the stability properties of the θ-method, consider a single variable

model problem:

du

+ »u = f , (15.13)

dt

with » > 0. This differential equation has the property that any perturbation to the

solution (for example induced by a perturbation of the initial value of u) decays

ˆ

exponentially as a function of time. Assume that u satis¬es the differential Eq.

(15.13) exactly, and let u be a perturbation of u, hence u = u + u. Consequently,

˜ ˆ ˆ˜

˜

the perturbation u must obey

˜

du

+ »˜ = 0.

u (15.14)

dt

If at t = 0, the perturbation equals u = u0 , the solution to this equation is

˜ ˜

u = e’»t u0 .

˜ ˜ (15.15)

Clearly, if » > 0, the perturbation decays exponentially as a function of time.

Application of the θ-scheme to the single variable model problem (15.13) yields

un+1 ’ un

+ θ»un+1 + (1 ’ θ) »un = θfn+1 + (1 ’ θ) fn . (15.16)

t

The one-dimensional convection-diffusion equation

268

˜ ˆ

Now, as before, if un is a perturbation of un , this perturbation satis¬es

u n + 1 ’ un

˜ ˜

+ θ»˜ n + 1 + (1 ’ θ) »˜ n = 0.

u u (15.17)

t

Clearly, the perturbation at t = tn + 1 can be expressed as

1’(1 ’ θ) » t

un + 1 =

˜ ˜

un . (15.18)

1 + θ» t

A

The factor A is called the ampli¬cation factor. To have a stable integration scheme

˜ ˜

the magnitude of un + 1 should be smaller than the magnitude of un , i.e. the pertur-

bation should not grow as time proceeds. Hence, stability requires |˜ n + 1 | ¤ |˜ n |,

u u

which holds if the ampli¬cation factor |A| ¤ 1.

Fig. 15.2 shows the ampli¬cation factor A as a function of » t with θ as a

parameter. For 0 ¤ θ < 0.5 the integration scheme is conditionally stable,

meaning that the time step t has to be chosen suf¬ciently small related to ».

In the multi-variable case the above corresponds to the requirement that, in case

0 ¤ θ < 0.5, » t should be small compared to the eigenvalues of the matrix A.

For 0.5 ¤ θ ¤ 1 the scheme is unconditionally stable, hence for any choice of

t a stable integration process results.

1.5

1

0.5

θ = 0.75

A 0

θ = 0.5

’0.5

θ = 0.25

’1

θ=0

’1.5

0 1 2 3 4 5

» ”t

Figure 15.2

Ampli¬cation factor A as a function of » t for various values of θ .

269 15.4 Spatial discretization

15.4 Spatial discretization

Following a similar derivation as in the previous chapter, the weak form is

obtained by multiplication of Eq. (15.1) with a suitable weighting function w,

performing an integration over = [ a, b], followed by an integration by parts:

‚u ‚u dw ‚u

dx + dx + dx = B,

w wv c (15.19)

‚t ‚x dx ‚x

where the right-hand side term B results from the integration by parts:

‚u ‚u

B = w( b) c ’ w( a) c . (15.20)

‚x ‚x

x=b x=a

Notice that no partial integration of the convective term has been performed. The

discrete set of equations, according to Eq. (15.19), is derived by subdivision of the

domain in elements and by discretization at element level according to

uh ( x, t) | = N T ( x) ue ( t) , wh ( x) | = N T ( x) we ( t) . (15.21)

∼

∼ ∼ ∼

e e

Note that the shape functions N are a function of the spatial coordinate x only and

∼

not of the time t. The nodal values of uh , at element level stored in the column

ue , however, do depend on the time t. Substitution of Eq. (15.21) into Eq. (15.19)

∼

yields

Nel

dN T

du

N N dx ∼e + wT

wT T

N a ∼ dx ue

∼e ∼e ∼

∼∼ ∼

dt dx

e e

e=1

dN dN T

+ = B.

wT ∼

c ∼ dx ue (15.22)

∼e ∼

dx dx

e

With

Me = N N T dx, (15.23)

∼∼

e

and

dN T dN dN T

Ke = N a ∼ dx + ∼

c ∼ dx, (15.24)

∼

dx dx dx

e e

Eq. (15.22) can be written as

Nel

due

+ K e ue = B.

wT M e ∼

(15.25)

∼e ∼

dt

e=1

After the usual assembly process this is written in global quantities:

du

+ K u = wT f ,

wT M ∼

(15.26)

∼

∼ ∼∼

dt

The one-dimensional convection-diffusion equation

270

where f results from the contribution of B (see Section 14.5). This equation has to

∼

be satis¬ed for all w, hence

∼

du

+ K u = f.

∼

M (15.27)

∼ ∼

dt

This is a set of ¬rst-order differential equations having a similar structure as Eq.

(15.11). Therefore, application of the θ-scheme for temporal discretization yields

1 1

M’(1 ’ θ) K un + f θ ,

M + θK un + 1 = (15.28)

∼ ∼ ∼

t t

with

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

∼ ∼ ∼

Clearly, in the steady case the set of equations, Eqs. (15.27), reduces to

K u = f. (15.30)

∼ ∼

Example 15.2 Consider the steady convection-diffusion problem, according to

du d du

=

v c ,

dx dx dx

with the following parameter setting:

= [ 0 1]

: x = 0 and x = 1

u

=…