<< ѕредыдуща€ стр. 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 стр.)ќ√Ћј¬Ћ≈Ќ»≈ —ледующа€ >>