Figure I. The spreadsheet usad to calculate CJRand h(CJR)
spreadsheets having a copying command and the ability to designate specific and relative cell addresses, calculations of this type can he accomplished with a few keystrokes. The column of cells E l . . .El4 contains the differences between D l . . Dl4 and D4. The column of cells F I . . F14 contains the squares of the differences in E l . . E14, and F4 contains the sauare root of the averaee of the souares of the differences. when this quality is rou;ded to the'proper number of signiiicant fieures. .. . in this case one. the value of AiC,IR) is 0.5 Application of eq 2 to this calculation leads toan involved aleebraic exnression (see Fie. 21. One can see that the confidence limitcalculated using eq 2 is identical to the direct spreadsheet method. Alternatlve Error Calculations Using Spreadsheet Programs One may proceed directly from eq 2 to a working equation for propagation of error that can also he automated on a personal computer spreadsheet. As suggested by a reviewer, onemay substitute for the partial derivative in eq 2, the close approximation AWAx for aFlax, etc. The most useful choice for Ax is A(x), a substitution that leads directly to eq 3
Implementation of the Direct Propagatlon of Error
Procedure on a Spreadsheet Program A calculation suitahle for the demonstration of the spreadsheet error propagation procedure is found in the method due to Clement and Dhsormes for the calculation of . ~ the heat capacthe heat capacity ratio for an i d e a l g a ~Here ity ratio of an ideal gas is determined from the manometrically determined values of PI, the gas pressure a t temperature TI; Pz, the gas pressure immediately after a reversible adiabatic expansion when the temperature is Tz; and Pa, the gas pressure after the gas has warmed back up t o temperature TI. The equation used to calculate the heat capacity ratio is
where
Reasonable confidence limits for all the independent variables are taken to be A(P1) = 0.2 torr, A(P2) = 0.1 torr, A(P3) = 0.2 torr. For purposes of simplicity and illustration, only the values of CJR and A(C,/R) will be calculated in this paper. The spreadsheet used tocalculate C,/R and X(C,/R) is shown in Figure 1. Cells A4, B4, and C4 contain PI, Pz, and P3,respectively. CellsA5, B5, and C5 contain the confidence limits for PI, Pz, and P3, respectively. The block of cells A7 . . C14 contains all the confidence limit comhinations. Since there are three independent variahles, there are eight possible comhinations of confidence limits. The calculated result CJR is in cell D4. The column of cells D l . . Dl4 are values of CJR calculated using the confidence limit combinations in the block A7 . . C14. For example, the value of C,/ R in cell D8 is calculated from the values in cells A8, B8, and C8; the value of CJR in cell D9 is calculated from the values in A9, B9, and C9; etc. Because most commonly available Rayleigh. Lord. The Theory of Sound, 2nd ed.; Dover: New York, 1945; Vol. 11, pp 15-23. 868
Journal of Chemical Education
The equation is easily automated on a spreadsheet program with the advantage that only 2n rather than 2" calculations are required. Moreover, eq 3 gives results identical to eq 1 when A(x) is small. An interesting exercise for the student is to allow A(x), Ah), . . . to increase and test how large the errors in the independent variahles can he before a significant divergence between the results of eq 1 and eq 3 are observed. Advantages ot the Dlreci Spreadsheet Method Slmp~icny Students can gain enough proficiency in the use of spreadsheets to perform this analysis in two lahoratory periods. The use of calculus is not required toestimate error since the same calculation that is used to obtain the average value of the dependent variable is used to obtain the confidence limit. Capability and Power To Explore the Significance of the Error Estimation Procedure After the spreadsheet is constructed, students can change the value of the confidence limit of any or all of the independent variahles and, after automatic recalculation of the spreadsheet, immediately see the result.
A Simple Numerical Method for Solving Complex Equilibria Jerald A. Devore California State University, Long Beach, CA 90840 A number of papers discussing methods of determining equilihrium concentrations in complex systems have appeared in this Journal (1-11). This simplex method (12,13) wherein the Gihhs free energy is minimized has become the standard technique for performing thermodynamic computations on reactive flow in nozzles. To a considerable degree, however, these methods are beyond the abilities of most undergraduates. In order to make the modeling of complex equilihrium systems accessihle to students of physical chem-
istry and scientific programming, we have adopted the simole iterative techniaue -Dro~osedbv Goldwasser (14) in khich initial concentrations df speoieiare repeatedly iefined so as to satisfy all equilibrium eauations. As an illustrative example, consider the phosphoric acid-phosphate system consisting of seven coupled equilibria:
H2PO;
PO:-
+ H,O e H8P04+ OH-
KJK,
+ H,O
KJK,
H,O
FA
HP0:-
+ OH-
* H+ + OH-
Kw
While the hydrolysis reactions can be written in terms of the acid dissociation equations and K,, their inclusion is necessary to achieve reasonable convergence rates in the pH regions where buffering occurs. The initial concentrations of various species in the bulk mixture are specified. All other concentrations are set equal to zero initially, except for [Hf] since no other sources are and [OH-], which are set to present. If [H3PO4Ii.[H2PO;li etc., are the molar concentrations a t the beginning of the ith iteration, the next iterate values computed from the Kt equilibrium equation are obtained by writing [ H , P O , ~= [H,POJ
-A
[ H , P O ; ~ = [H,PO;)~
+A
[H+r+l= [H+p+ A
The adjustments to concentrations are then determined by substituting the above relations into the equation for KI.
Terms containing powers of A greater than 1 are dropped giving
SUBROUTINE CALCICI Subroutine CALC iteratively determinsr concsntrationr tor the p b ~ p h a l e~ystem,given a Eet of initial candtione. C(1)=[H3P04]. C(2) = [H2P04(-)I C(3) = IPHOR(2~)) C(4) = IP04(3~11,C(5) = [H(+)]. C@I=[OH(-11 DOUBLE PRECISION C(6),CTEMP(6).K(3).KW,DELTA DOUBLE PRECISION RETT RESID. SET
.
.
~~=1.0$-14 The parameter SET detsrmlnes the convergence criteria. SET=t.ODB RESID=t.O DO WHiLE(RES1DGTSET) CTEMP values store mncsntralion from the previous iteration. DO 10 1-1.6 CTEMP(I,=C(I) The three dlsso~lallonreanions are solved in turn. DO 20 J=1.3 DELTA-(K(J)+C.Jl- C(J + 1)*Cl51)IlC(J+ 11+ Cl51 t If any mncentratron 8s negat8ve. skip to next squatton. IF(C(J) - DELTA.LE.O.OR.CtJ+ 1) + DELTALEO ORC(5I + DELTA.LE.O) GOT0 20 C(J)=C(J)- DELTA C(J + 1 I = C(J + I ) + DELTA C(5) = C(5) + DELTA CONTINUE The three hydrolysis reanlons are solved in turn. DO30J-1.3 DELTA=(KWIK(J),C(J + 1) - C(J).C(6))I(C(J) +
Water dissaciaton equation is solved DELTA = (KW - C(5).C(6))lC(5)+ C(6)). It [Hi+)] 01 [OHi~)]is negative, sk~pto next iteration. IF(C(5)+ DELTALEOORC(6) + DELTALEO) GOTO40 C(5)= C(51r DELTA C(6) = C(6) + DELTA CONTINUE The relative change in concentration of each specie* from it$ value at the previous iteration is caicuiated and set sqvai to RETT The largest of these 1s equafed with RESID whlch iscornpared with the conueigence criterion at each hation. Concentrat8ons are also tested lor zero values as a precaution against possible program inlerupllan. IF(CI1)NE.O) RESiD-ABS(CTEMP(1)- C(I))IC(l) DO 50 1b2.6 IF(C(1)NE.O)M E N RETT = ABS(CTEMP(I1- C(I))IC(Il IF(RETT.GT.RESID)RESID = RETT END IF CONTINUE END DO RETURN END
ComposMon and pH Valves at Selected Tltratlon Polnts Figure 1. Subroutine for iteratively solving phosphate equilibrium equations. OH- (mL)
[HsPO+]
0.00 15.15 20.10 26.20 39.75 46.00
0.7655-01 0.1319e-01 0.746%-04 0.3558-06 0.6195-10 0.1467e15
[ H s P O ~ ~ ] [HP0,2-] 0.2335-01 0.4907-01 0.5500e-01 0.2773~01 0.5144e-03 0.6666-06
0.6340e-07 0.1626e-05 0.3614e03 0.1927%01 0.3806e01 0.2700e-01
[PO?] 0.1141e-17 0.357%-15 0.1573-10 0.8869-07 0.1867-04 0.7245e-02
pH 1.63 2.72 5.02 7.04 9.07 11.81
A=
KL[H3PO4Ii - [HzPO;li[Htli [H,PO;I'+ [Hip
+X,-
Each equilibrium equation is treated similarly in turn so that the concentrations of most species are adjusted several times during one iterative cycle &er the entire set of equations. If a particular reaction yields a negative concentration. that reaction is s k i.o.~ e duntil the next iteration. All reactions will eventually give positive concentrations as eauilibrium is ao~roachedand the deltas become sufficientl ~ s m a l lMass . baiance is maintained automatically with this procedure if no separate phases are present. A FORTRAN 77 subroutine for solving the phosphate eauations is shown in Fieure 1. In this examole. convergence is'achieved when therelGtuechangeof all species is leisthan 10.' between iterations. The"housekeeoinr" shell Droaram. which is capable of performing automaiic &rations, or simDIY c o m ~ u t i n ethe eauilihrium values for a s~ecifiedhulk mixture and optionaliy adding specified amounts of acid or base. has been omitted for brevitv. The table shows several selected points from a simulatedtitration of 0.1 M H3P04 with 0.125 M OH-. Concentrations of all species are dis-
-
Figure 2. Simulated titration of 25 mL of 0.1 M H3P04wlth 0.125 M OH-
Volume 65 Number 10 October 1988
869
played allowing the student to appreciate immediately the complex buffering that occurs a t higher pH values. Figure 2 illustrates the pH titration curve as a function of the volume of base added. In this example, the program was set to print out pH values a t about 0.15 intervals. Base was added in 0.15-mL increments and the equilibrium composition determined after eachaddition. The iteration count also shown in Figure 2 represents the total number required from the last point printed, which in some regions requires several additions of base. Only seven iterations are required to equilibrate the original phosphoric acid solution, whereas 410 iterations are required to determine the point a t 25.05 mL of base delivered from the previous point a t 23.85 mL. By contrast, 93 iterations are required if 25.05 mL of base is added immediately to the original equilibrated solution. During a titration simulation, the number of iterations required depends both on how many OH- increments are required and on how rapidly the composition is changing from point to point. One advantage of this method is that convergence is achieved regardless of the order in which the equilibrium equations are solved. Roughly a 5% difference in the number of iterations required for convergence is observed if the equations are solved in the reverse order from that given in Figure 1.During a titration this difference can accumulate if there are several additions of base between listed points. At low pH values, the order given in Figure 1is more efficient, whereas at high pH values the reverse order requires fewer iterations. Double precision arithmetic is recommended for subroutine CALC when operated on machines with less than 60-bit word lengths. Figure 3 illustrates concentrations of the four phosphate species over the entire titration. The sharply peaked H2PO; and HPOq- concentration curves are interesting for the student to observe as this behavior is only occasionally mentioned in discussions concerning titrations. Figure 4 shows the distribution diagram for phosphate species as a function of pH and is a useful tool for preparing buffered solutions. The contents of this diagram were generated by making several small changes in the automatic titration routine. A more ambitious group project using this technique was successfully undertaken by chemical engineering students studying physical chemistry. The performance of gasohol compared to gasoline was studied using a reversible Otto thermodynamic cycle. This type of simulation is quite similar to calculating theoretical performance parameters of rocket engines and requires determining flame temperatures and shifting equilibrium under various thermodynamic constraints. A minimum of about 25 equilibrium equations are required for a combustion problem consisting of the elements C, H, 0 , and N in order to account for high temperature species. The elements present as reactants are assigned to the principal combustion products H20, CO, C02, and C(s) initially to start the calculation. The amount of solid carbon present is tracked using a separate carbon mass balance equation. The simplicity of this iterative method allows the student to comprehend the logic and create an application without great difficulty. Convergence to the equilibrium composition is virtually assured given sufficient equilibrium reac-
870
Journal of Chemical Education
Figure% Composition of phosphatespeciesduring simulated titration of 25 mL of 0.1 M H3P0, with 0.125 M OH-
Figure 4. Disbibution of specles diagram f w phosphoric acid.
tions to describe the chemistry of the system. The order in which equilibrium equations are solved is important only in that a judicious choice may speed convergence somewhat. Separate mass balance equations are needed only if the system contains more than one phase. 1. Bard,A.J.;King,D. M. J . Cham.Educ. 1965.42.127. 2. Emery, A. R.J. Chem. Edue. 1965,42,131. 3. Zajicek, 0. T.J Chem. Educ. 1965,42,622. 4 . Eberhart,J.G.J.Chem.Edue. 1965,42,624. 5. Senesf, D.J. Chsm. Educ. 196S,42,625.
6. Bard,A. J. J. Chern. E d u . 1965.42.626. 7. Stone. E. E. J . ChemEdue. 1966.43.241. 8. Haglund, E.: Moan, D.: Flynn, J. J . Chem. Educ. I966,43,582. 9. Shonfeld, E. J . Chrm.Edue. 1968.46, 173. 10. Emersan,J.H.J C h e m E d u c . 1978,55,641. 11. Feens0s.T. P. J. Chem.Educ. 1979,56,104. 12. White, W . B.:Johnson. S. M.: Dsnfzig, G.B. J.Chem. Phya. 1958 13. Boynton.F.P. J. Chem Phys 1960.32.1380. 14. Go1dwasser.S. R.Ind.Ena. Chem. 1959,51,595.
