<<

. 54
( 67 .)



>>

0 0.2 0.4 0.6 0.8 1
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
=…

<<

. 54
( 67 .)



>>