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.