Environmental Engineering Reference
In-Depth Information
Imagine that we don't know the analytical solution, or would not have the
possibility to evaluate it; maybe we have a supermarket calculator only without
access to the exponential function. Are there other means to obtain a solution value
for a given time
t
?
Luckily there are numerical methods. For the current problem we will outline the
simplest numerical method. The idea is to approximate the differential term on the
left side of the (
21.1
) by a finite difference term. The numerical method is thus
termed '
finite differences
' (FD). Here the finite difference looks like this:
@c
@t
c
ð
t
þ
D
t
Þ
c
ð
t
Þ
D
t
(21.4)
D
t
is the so called timestep, and we can use the approximation for any
t
that we
like. Geometrically the approximation can be visualized by taking the secant
through the two points (
t
,c(
t
)) and (
t
+
D
t,c
(
t
+
D
t
)) instead of the tangent at the
point (
t
,c(
t
)). From (
21.1
) we obtain:
c
ð
t
þ
D
t
Þ
c
ð
t
Þ
D
t
¼
l
cðtÞ
(21.5)
Assume that
c
(
t
) is known. The we have to resolve equation for the unknown
c
(
t
+
D
t
), and obtain:
cðt þ
D
tÞ¼ð
1
l
D
tÞcðtÞ
(21.6)
Thus we see that we have to multiply the concentration at time
t
with the
factor 1
l
D
t
, to obtain an approximation of the concentration at time
t
+
D
t
.
The factor depends, aside from the given material property
l
, only on the numerical
parameter
D
t
.
We utilize the formula (
21.6
) in the following 'timestepping' approach. We start
at initial time
t ¼
0, at which the concentration is known (formula (
21.2
)) and apply
formula (
21.6
). In that way we obtain an approximate value for the concentration at
time
D
t
:
cð
D
tÞ¼ð
1
l
D
tÞc
0
(21.7)
Now that we have a value for
c
at
t ¼
D
t
, we can apply the formula (
21.6
) once
more to obtain a value for
c
at time
t ¼
2
D
t
:
cð
2
D
tÞ¼cð
D
t þ
D
tÞ¼ð
1
l
D
tÞcð
D
tÞ
(21.8)
In that manner we use the formula (
21.6
) over and over again until we cover the
entire time interval in we are interested and reach the final time instant that we want.
The entire procedure is implemented in the following M-file. Starting time is
zero and final time is given by the variable
Tmax
. After the initialization of all
