APPLICATIONS OF TECHNOLOGY IN TEACHING CHEMISTRY
An On-Line Computer Conference
June 14 TO August 20, 1993
PAPER 13
Finite Difference Solution of the Diffusion Equation on a Spreadsheet
Douglas A. Coe, Montana College of Mineral Science and Technology,
Butte, MT 59701 DACOE%MTVMS2.MTECH.EDU
SCHEDULE: Short questions on this paper: July 28, 1993
Discussion of this paper: Aug. 6 through Aug. 9, 1993
I ABSTRACT
Physical chemistry topics involving second order partial
differential equations often receive only cursory treatment in
undergraduate curricula because analytic solutions of these equations
are readily obtainable only for the simplest of boundary conditions.
This paper will describe the solution of the diffusion equation on a
spreadsheet, using the method of finite differences. The solution will
be illustrated for the diffusion of Ar at 25 C through a 0.100 mm thick
polyvinyl acetate membrane. Graphical comparison of the spreadsheet
and analytic solutions show them to be virtually identical.
II PAPER
Paragraph A:
This paper illustrates how partial differential equations
that govern flow phenomena, specifically the diffusion equation,
can be solved on a spreadsheet using the method of finite
differences. The diffusion equation in one dimension, when the
diffusion coefficient is not a function of distance diffused or
time, is (1):
dC/dt = D d^2C/dx^2 1
In this equation C is the concentration of the diffusing species, t is
the time the species has diffused, x is the distance the species has
diffused, and D is the diffusion coefficient. dC/dt is the first partial
derivative of concentration with respect to time. d^2C/dx^2 is the
second partial derivative of concentration with respect to distance.
Paragraph B:
In the finite difference method the time and distance are not
allowed to vary continuously, but rather are incremented in small
discrete steps. The time and distance derivatives in eq 1 are
represented by Taylor series expansions ignoring, respectively, terms
higher than linear and quadratic in each expansion. This leads to the
finite difference approximation of the diffusion equation:
C(n+1,m) = C(n,m)+(D*Del(t)*(C(n,m+1)-2*C(n,m)+C(n,m-1))/(Del(x))^2 2
C(n, m) is the concentration of the diffusing species at time step n and
distance step m. Del(t) is the size of the time step delta t, while Del(x)
is the size of the distance step delta x. Equation 2 allows the
concentration at the new time step, n+1, to be calculated from
concentrations at the previous time steps, n.
Paragraph C:
The diffusion of Ar in the gas phase at 25øC through a 0.100 cm think
polyvinyl acetate membrane (2) will be used to provide a specific
example of the finite difference solution of the diffusion equation on a
spreadsheet. A portion of the spreadsheet used in this example is
shown in Figure 1. The total membrane thickness of 0.100 cm is divided
into 20 distance steps of 0.0050 cm each, incremented horizontally in
the spreadsheet in cells B16 to V16. Time is incremented vertically in
the spreadsheet in steps of 1000 sec (0.28 hrs) in cells A18 to A258.
The initial condition that the membrane contain no Ar is imposed in the
spreadsheet by requiring that all 20 cells in the membrane at time zero
show an Ar pressure of zero bar, as indicates by the values of 0.00 in
cells C18 to V18. The boundary conditions of a constant pressure of
0.40 bar at the left membrane surface and a vacuum at the right
membrane surface are imposed in the spreadsheet by requiring that
cells B18 to B258 contain the value 0.40 and cells W18 to W258 contain
the value 0.00. The diffusing Ar pressures in the membrane are given
by values in the cell range C18 to V258. As a typical example of the
formulas used to calculate the Ar pressures in the membrane (see eq 2),
consider the formula in cell D25:
IF($E$13=1, 0, D24 + $E$11 * (E24 - 2*D24 + C24)) 3
The IF statement allows the membrane to be initialized with zero Ar
pressures by calculating the spreadsheet when cell E13 contains a value
of 1. Replacing the 1 in cell E13 with any other value, e.g. 0, and
recalculating the spreadsheet, then allows the diffusion calculations to
be carried out from this initial state. The value of the product
D*Del(t)/(Del(x))^2, on which eqs 2 and 3 depend, is explicitly included in
cell E11, since the finite difference solution to the diffusion equation
exhibits unstable oscillations, when this product has a value greater
than or equal to 1/2.
Paragraph D:
The calculations were continued until the steady state pressure
profile of Ar in the membrane is no longer changing with time, at which
point the Ar pressure in the membrane should vary linearly with
distance diffused. Figure 2 shows that a steady state pressure profile
was not reached until the Ar had diffused for 66.7 hours. The
restriction that D*Del(t)/(Del(x))^2 < 1/2, with a diffusion coefficient of
1.08 x 10^(-8) cm^2/sec and a Del(x) of 0.0050 cm requires Del(t) to be
less than 1157 sec. The calculation illustrated in Figure 2 uses a Del(t)
of 1000 sec and, accordingly, needs 239 time steps to reach steady
state. Since each time step corresponds to a row in the spreadsheet,
the memory requirements of the spreadsheet are quite large.
Paragraph E:
The analytic solution to the diffusion equation, where the pressure
of the diffusing species is initially zero in the membrane and is also
held at zero at the right membrane surface is (3):
P = Pl + (Pr-Pl)*x/L + (2/pi)*SUM(n,(1/n)
*(Pr*COS(n*pi)-Pl)*SIN(n*pi*x/L)*e^(-D*n^2*pi^2*t/L^2)) 4
In this equation P is the pressure of the Ar within the membrane, Pl is
the pressure of the diffusing species at the left surface of the
membrane, Pr is the pressure of the diffusing species at the right face
of the membrane, L is the membrane thickness, e is the base of the
natural logarithms ,and the other quantities are the same as previously
defined. The notation SUM(n,argument) means to sum the argument over
the summation index, n. Figure 3 compares the spreadsheet and analytic
solutions of the diffusion equation with solid curves representing the
analytic solutions and data symbols the spreadsheet solutions. At 4.44
hours, early in the diffusion process and well before steady state has
been reached, there is no discernible difference between the analytic
and spreadsheet solutions. At 66.7 hours, when steady state has been
effectively reached, the agreement between the analytic and
spreadsheet solutions is still in general very good, but does differ
slightly near the right surface of the membrane. This difference is
attributable to neglecting higher order terms in the Taylor series
expansions used to derive eq 2.
Paragraph F:
The calculations were done on a 486 microcomputer with a math co-
processor using Borland's Quattro Pro 123 spreadsheet (version 4.0).
The diffusion spreadsheet uses 342 kilobytes of RAM. The diffusion
spreadsheet takes ~20 seconds to load into memory from disk, ~20
seconds to save to disk, and ~10 seconds to recalculate from the initial
condition in which the pressure of Ar is zero in the membrane.
Paragraph G:
Construction of this spreadsheet will illustrate the finite
difference solution of partial differential equations and the inherent
tradeoff between calculational accuracy and the speed and memory
requirements in numerical calculations. Once constructed the
spreadsheet can be used by the student to develop an intuitive feeling
for diffusion by studying the effects on diffusion of varying the
diffusion coefficient and concentrations at the membrane surface and
within the membrane on the time and position development of diffusion
profiles.
III REFERENCES
1. Crank, J. The Mathematics of Diffusion; Clarendon Press: London,
1956; p. 4.
2. Meares, P. J. Am. Chem. Soc. 1954, 76, 3415-3422.
3. Crank, J. The Mathematics of Diffusion; Clarendon Press: London,
1956; p. 47.
IV QUESTIONS
1. Are these sort of exercises of any pedagogical value?
2. What is the difference between what the student learns if they
have to construct the spreadsheet versus being given a working
spreadsheet of this model?
3. Does the typical undergraduate chemistry student have enough
knowledge of spreadsheets to build this model. Is this class or
institution dependent? Is exposure to second order partial
differential equations a prerequisite?
4. What is the relative educational value of exposing students to (1)
the diffusion equation, (2) numerical solutions of differential
equations, (3) an advanced spreadsheet exercise?
5. Some effort is required by the student to construct the
spreadsheet described in this paper. Is exposing students to
numerical solutions of the diffusion equation on a spreadsheet
worth the effort?
V FIGURES
Figure 1. This is a portion of the spreadsheet used in the finite
difference solution of the diffusion equation. The file
FIGURE1.TXT was originally created in Wordperfect 5.1 to look
like the spreadsheet and then saved as an ASCII file.
FIGURE1.TXT is 2 KB.
Figure 2. This figure gives the spatial distribution of pressure of Ar
diffusing through the membrane at different times. The file
FIGURE2.PIC was originally saved as a .PIC file from the
spreadsheet Quattro Pro 123. The .PIC file was then converted
to a .PCX file, using PCXDUMP.EXE. The .PCX file was converted
to a 640x480 .GIF file using PICEM.EXE. Finally the .GIF file was
encoded for network transfer as an .UUE file, using
UUENCODE.EXE. FIGURE2.GIF is 9 KB.
Figure 3. Comparison of the finite difference (solid lines) and analytic
solutions of the diffusion equation at 4.44 hours (squares) and
66.67 hours (asterisks). The file FIGURE3.PIC was originally
saved as a .PIC file from the spreadsheet Quattro Pro 123. The
.PIC file was then converted to a .PCX file, using PCXDUMP.EXE.
The .PCX file was converted to a 640x480 .GIF file using
PICEM.EXE. Finally the .GIF file was encoded for network
transfer as an .UUE file, using UUENCODE.EXE. FIGURE3.GIF is
9 KB.
VI BINARY FILES
A 400 KB binary file of the Quattro Pro 123 (version 4.0)
spreadsheet DIFFUSE.WK1 has been included with this paper and can be
obtained by anonymous ftp from:
info.umd.edu info/Teaching/ChemConference/Papers.