)V1ESI_J RETURNING MATERIALS: P1ace in book drop to LIBRARIES remove this checkout from J-uqu-u-L your record. FINES wil] be charged if book is returned after the date stamped below. APPLICATION 017 um GENERAL mnsron'r THEORY T0 NONISO'I'HERIIAL GALVANIC CELLS By Bruce Kennedy Borey, Jr. A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Chemistry 1984 @oopmcmm BRUCEKENNEDYEDREY, JR. 1934 ABSTRACT APPLICATION OF THE GENERAL TRANSPORT THEORY TO SOLID STATE GALVANIC CELLS By Bruce Kennedy Borey, Jr. The goal of the research described here is the solution of the set of partial differential equations from nonequilibrium thermodynamics, hydrodynamics, and electrostatics which describes ion transport through a variety of energy conversion devices. The solution is obtained both numerically and analytically under a variety of boundary conditions. This permits simulation of many different electrode[electrolyte/electrode systems. The electrolyte spans the solution space and is assumed to be a continuous phase, incompressible, free of chemical reactions, and subject only to conservative external forces. The solutions to the differential equations are in the forms of both numerical finite difference simulations and approximate analytical expressions. The importance of electrolytic properties in determining the cell efficiency is evaluated. Inclusion of temperature as a dependent variable allows examination of the feasibility of using temperature gradients to enhance economic performance. Contact vith experiment is made through the numerical calculation of the ohmic voltage drop across the electrolyte. Extrapolation of the curve of ohmic drop versus current yields a value for the bulk conductivity in good agreement with the value previously obtained from a.c. measurements. To gain information concerning overpotentials, a.c. circuit theory is employed. Certain circuit components. (resistors, capacitors). combine to form circuits whose behavior mimics that of real systems. Measurements of cell impedance: taken over a wide range of frequencies are analyzed in the complex plane. The analysis leads to estimates of the effect of overpotentials upon cell efficiency. Mass transport transient behavior is examined and a simple formula is derived for calculating the time required for establishment of the steady state. To Marilyn. Brucie, Sara, and Erica ii ACKNOWLEDGEMENTS I wish to thank the department of chemistry for the financial support I received from teaching assistantships. I also wish to thank Mobil Oil Corporation for awarding me a fellowship through the department of chemistry in 1982. I am grateful to Professor Frederick H. Horne for his encouragement and patience throughout my stay at Michigan State University. I appreciate his many helpful suggestions in the preparation of this dissertation. I am also grateful to the other members of my committee, Professor Katharine C. Hunt for acting as second reader and making many helpful suggestions, and Professors Thomas J. Pinnavaia and Harry A. Eick for agreeing to serve on my committee on such short notice. I would like to express my most sincere appreciation to Dr. Tom Atkinson for his assistance in utilizing the computer facilities throughout my stay at Hichigan State University. I am especially grateful for the many hours he spent helping me use the spinwriter which produced the final copy of this dissertation. I am also grateful to Mr. Dwight Lillie for his valuable assistance in providing me with the computing facilities necessary for the completion of this dissertation. I am indebted to other members of my research group, Dr. John 111 Lockey, Dr. Richard Rice, Mr. Tom Troyer. and Mr. Yuan In. for their helpfulness and friendship over the years. I would also like to thank the many students and faculty who have given me their friendship and from whom I have learned much. I thank my parents for their love and encouragement. and for providing me with such a great source of pride. I would also like to thank my father and mother-in-law, Ir. and Hrs. G. Lewis Potter for their kindness and generosity. but especially I thank them for their daughter. lost of all I thank my wife. larilyn, who endured years of hardship enabling me to complete this work. She never lost faith, and without her I could never have completed this work. iv TABLE OF CONTENTS Chapter Page LIST OF TABLES . . . . . . . . . . . . . . . . . . ix LIST OF FIGUMS. C O C I O O O O C O C O O O O O O x CHAPTER 1 - MACROSCOPIC ANALYSIS OF ENERGY CONVERSION DEVICES . . . . . . . . . . 1 A. Introduction . . . . . . . . . . . . . . . . 1 B. Objectives. . . . . . . . . . . . . . . . . 2 C. Plan of the Dissertation . . . . . . . . . . . . 3 CHAPTER 2 - ION TRANSPORT EQUATIONS . . . . . . . . . . . 6 A. Introduction . . . . . '. . . . . . . . . . . 6 B. Nonequilibrium Thermodynamic Equations . . . . . . . 6 C. Continuity Equations . . . . . . . . . . . . . 11 D. Electrostatic Equations . . . . . . . . . . . . 11 E. 'Combined System of Equations. . . . . . . . . . . 13 F. Solving the Equations . . . . . . . . . . . . . 15 CEAPTER 3 - BOUNDARY CONDITIONS . . . . . . . . . . . . 13 A. IntIOduction O O O O O C O O O O I O O C O I 18 Chapter CHAPTER Page Physical Interpretation of the Bound.” condi t ion' 0 O I O O O O O O O O O O O 21 Applications of the Boundary Conditions . . . . . . . 27 4 - STEADY STATE TRANSPORT IN OPEN SISTERS . . . . . . 34 Introduction . . . . . . . . . . . . . . . . 34 Description of the System. . . . . . . . . . . . 34 Steady State Equations. . . . . . . . . . . . . 37 Solution to the Isothermal Problem. . . . . . . . .. 39 Determination of Efficiency . . . . . . ‘. . . . . 47 Results and Discussion. . . . . . . . . . . . . 49 1. Introduction. . . . . . . . . . . . . . . 49 2. The Hixed Conductor System . . . . . . . . . . 50 3. Calculation of Efficiency . . . . . . . . . . 53 5 - STEADY STATE NONISOTHERIAL TRANSPORT IN 0m sysms O O I I O C O O O O I O O O 6 1 Introduction . . . . . . . . . . . . . . . . 61 Temperature Distribution . . . . . . . . . . . . 62 Solution to the Nonisothermal Problem. . . . . . . . 66 vi Chapter Page D. Efficicncy. O O O O O O O O O O O O O O O O 68 E. Results and Discussion. . . . . . . . . . . . . 73 CHAPTER 6 - NUNERICAL APPROACH. . . . . . . . . . . . . 79 A. Introduction . . . . . . . . . . . . . . . . 79 l B. Finite Difference Formulation . . . . . . . . . . 82 C. Numerical Results . . . . . . . . . . . . . . 85 CHAPTER 7 - THE TIRE DEPENDENT PROBLER . . . . . . . . . . 100 A. Introduction . . . . . . . . . . . . . . . . 100 B. Formulation of the Problem . . . . . O. . . . . . 101 C. Solution of the Soret Problem . . . . . . . . . . 106 CHAPTER 8 - FUTURE IORK . . . . . . . . . . . . . . . 115 APPENDIX A - TRANSPORT COEFFICIENTS . . . . . . . . . . . 118 APPENDIX B - FIRST ORDER SOLUTION FOR THE STEADY STATE ISOTHERNAL PROBLEM . . . . . . . . . . 121 APPENDIX C - DINENSIONLESS VARIABLES USED FOR COIPUTER SIMULATION . . . . . . . . . . . . 129 APPENDIX D - JOULE HEATING . . . . . . . . . . . . . . 131 TRANSPORT EQUATIONS IN FINITE Dlmcz F0“ 0 I O O O O O O O O O O O 133 APPENDIX E 'vii Chapter Page APPENDIX F - LISTING OF PROGRAN "IONFLO" . . . . . . . . . 137 mm 0 O O O O O O O O O O O O O O O O O O 1 80 viii Table LIST OF TABLES LEAST SQUARES VALUES OF THE EMPIRICAL CONSTANTS GIVEN IN THE EQUATIONS (“v-OI)=a+bp,n= c + dp, XND R = e + fp FOR VARIOUS VALUES OF THE CURRENT. . . . . . LEAST SQUARES VALUES OF TUE EMPIRICAL CONSTANTS GIVEN IN I] = a + blO‘AT/T] FOR VARIOUS VALUES OF THE CURRENT. . . . . . . . . . COMPARISON OF Ad OBTAINED FROM INTEGRATION OF EQ. (4-30) ANS THE COMPUTER SIMULATION RESULTS FOR THE DOPED SOLID ELECTROLYTE ceoz + C30 a o e O O O 0 0 0 COMPARISON OF THE VALUES OF R AND R ELECTROLYTE cw2 + CaO . . . . . COMPARISON OF POWER DENSITIES WITH AND WITHOUT OVERPOTENTIALS . . . . . . . . ix MEASURED AT HIGH TEMPERATURES FOR THE PED S LID Page 56 78 86 99 99 LIST OF FIGURES Figure 3-1 INTERFACIAL REGION AT A PHASE BOUNDARY . . . . 3-2 DEPENDENCE OF THE KINETIC RATE CONSTANTS ON a. 1. O O O O O O O O O O O O O 0 3-3 CIRCUIT REPRESENTATIONS 0F SOLID SOLID INTERFACES. (a) PARTIALLY BLOCKING ELECTRODE. (b) IDEALLY POLARIZED ELECTRODE. (c) IDEALLY DEPOLARIZED ELECTRODE. . . . . . 3-4 SUBDIVISION OF THE ELECTROLYTE INTO THREE PSEUDOPHASES WITH PSEUDOPHASE BOUNDARIES SHWN ATX =iL/2. o o o e o e e o o o 4-1 NONISOTHERMAL GALVANIC CELL. . . . . . . . 4-2 DEPENDENCE OF EFFICIENCY ON CURRENT DENSITY. (a) OHM'S LAW BEHAVIOR. (b) BEHAVIOR PREDICTED BY EQ. (4-57) 0 O O O O O O O O O O 0 4-3 DEPENDENCE OF R 0N CURRENT DENSITY . . . . . 5-1 DEPENDENCE OF THE EFFICIENCY ON CURRENT DENSITY. (a) O'AT < O, (b) O‘AT = o. (c) O‘AT > O . . . 6-1 NONUNIFORM SPACE AND TIME GRID . . . . . . 6-2 DEPENDENCE OF THE OHNIC VOLTAGE DROP 0N CURRENT DENSITY. . . . . . . . . . . . 6-3 EQUIVALENT CIRCUITS AND THEIR CORRESPONDING ADMITTANCE PLOTS. THE CONDUCTANCE (G) IS PLOTTED VERSUS THE.SUSCEPTIBILITY (B) . . . . 6-4 WARBURG IMPEDANCE PLOT . . . . . . . . . Page .20 24 29 31 36’ 58 60 76 81 89 92 95 Figure ' Page 6-5 DEPENDENCE OF THE EFFICIENCY ON THE CURRENT DENSITY. (+) NEGLECTING OVERPOTENTIALS. (o) INCLUDING W8 O C I O O O O O O O I O O O O O 98 xi CHAPTER 1 MACROSCOPIC ANALYSIS OF ENERGY CONVERSION DEVICES A. Introduction A fuel cell is a device for converting chemical energy (fuel) to electrical energy without the use of moving mechanical parts. Extensive studies on a wide variety of fuel cells using different electrodes, electrolytes. and geometries have been carried out (Etsell and Flengas [1970] and Foley [1969]). There are generally three sources of power loss associated with the steady state Operation of such devices, and they are the central foci of research in the field of energy conversion today. These three sources of power loss are classified as activation, mass transfer. and ohmic overpotentials. Of these the first two arise as a consequence of the kinetic limitations of the electrolyte/electrode interface involved. The ohmic overpotential is a result of ionic transport across the bulk or electrOlyte phase. A theoretical treatment of these phenomena is possible using the nonequilibrium .thermodynamic formalism. We are interested in the performance characteristics of the fuel cell in the steady state where the production of entrOpy has reached a minimum (Onsager [1931]). The combined nonequilibrium. hydrodynamic, and electrostatic set of equations can be used for many purposes. Rice [1982] uses these equations to describe the steady state isothermal transport] of ions through a membrane. Tasaka. Morita. and Nagasawa [1965] have investigated membrane potentials in nonisothermal systems. Several authors have studied the transport of vacancies, electrons, and holes in solid state electrolytes using the macroscopic formalism (Howard and Lidiard [1963]. Ohachi and Taniguchi [1977]. Cheung, Steele, and Dudley [1979]. Dudley, Cheung, and Steele [1980], and Weppner and Huggins [1977]). This dissertation deals with the theoretical study of transport phenomena in the steady state Operation of fuel cells, and how independent variables such as temperature, current, and external load can be adjusted to optimize performance. B. Objectives The principal objective of this work is to relate in a general way the steady state cell performance to electrolytic properties and externally adjustable parameters such as load voltage and temperature. This relationship is derived from the macroscopic equations of nonequilibrium thermodynamics. hydrodynamics. and electrostatics which govern the transport of matter. and is independent of the events occuring at the electrode electrolyte interface. This phase of the investigation thus deals only with ways of minimizing the ohmic overpotential losses suffered by any operating cell. Inclusion of a temperature gradient in the formulation of the problem is undertaken in order to assess the feasibility of using the nonisothermal properties of an electrolyte for the purpose of increasing power output. The assumption of electroneutrality or the neglect of regions of space charge in many of-the treatments cited above seems rather unrealistic in light of the fact that one expects a buildup of charge in the region of electrolyte near the charged surfaces of the electrodes. Provisions are therefore made for its inclusion in the formulation of the boundary conditions of the problem. and the importance of non-electroneutrality in electrolytic behavior is evaluated. The evaluation of nonohmic overpotential losses is undertaken with the aid of equivalent circuit anaIOgies and simple a.c. circuit theory. The goal here is to assess the relative importance of these sources of power loss. Power losses due to concentration and activation overpotentials are dependent upon the kinetic parameters of the electrode/electrolyte interface and are related to kinetic rate constants derivable from the macroscOpic transport equations (Horne, Leckey. and Perram [1981]). C. Plan of the Dissertation In Chapter 2 we begin with a review of the fundamental equations of nonequilibrium thermodynamics. hydrodynamics. and electrostatics used to formulate the macroscOpic transport equations. The last part of Chapter 2 deals with the separation of the starting equations into a part due to bulk behavior and a part due to nonelectroneutrality. and the specialization of these equations to the steady state. Chapter 3 is devoted to a discussion of the boundary conditions and in particular their application to Open systems in which the irreversible transfer of charge is taking place. In Chapters 4 and 5 analytical solutions to the steady state equations are derived for the transport of ions in galvanic cells under a variety of boundary conditions. These solutions highlight the separation of the bulk behavior of the system from behavior due to non-electroneutrality. Also the role of the nonisothermal properties of the electrolyte in enhancing cell performance is discussed (Borey and Horne [1982]). The importance of higher order terms in the. analytical solutions is considered and the performance of the cell as a function of current is examined. Chapter 6 begins with a brief discussion of the numerical formulation of the transport equations (Leckey [1981]) and their modification and adaptation to include nonisothermal terms. A comparison of numerical results with analytical results is made in an attempt to verify the accuracy of the latter, and to make contact with experiment. Simple a.c. circuit theory in conjunction with equivalent circuit representations is invoked to extract information about interfacial parameters such as contact resistances. This information is then ‘Used to estimate the relative importance of nonohmic overpotentials. In Chapter 7 we briefly examine the time dependent problem. The accuracy of the solution is verified by comparison to the steady state solution. and the transient behavior of the system is examined. The transient behavior of the system is taken to be the behavior over the period from the closing of the circuit to the establishment of the steady state. A simple formula is derived from which the time required to complete this transient phase can be estimated. This transient behavior in many applications may be of great importance since operating periods may be of comparable or shorter length than this time range, making steady state operation the exception rather than the rule (Jasinski [1967]). CHAPTER 2 ION TRANSPORT EQUATIONS A. Introduction The macroscopic equations used of the fields hydrodynamics. thermodynamics. Details of the phenomenological equations of to describe ion transport come from electrostatics. and nonequilibrium assumptions which lead to the nonequilibrium thermodynamics are available from many sources (deGroot and Mazur [1962], Haase [1969]. and Horne [1966]). The systems dealt with in this dissertation are one-dimensional, binary electrolyte-solvent mixtures in which pressure effects are negligible. B. In a charged system there (solvent, cation, and anion) fluxes and four diffusion the chosen as independent diffusion fluxes defined by 11 3 ci(vi-vo), 131,2. coefficients. Nonequilibrium Thermodynamic Equations are three independent velocities and therefore two independent diffusion With the cation and anion species (1 and 2 respectively) the Hittorf (2-1) where V1 is the velocity of species i and Va is the solvent velocity. are. according to Onsager [1931]. linear combinations of the gradients in electrochemical potential ii and temperature T. -Ji = 2 1ij[afij/ax]T + liq[alnT/axJ. (2-2) Here x represents the space coordinate and the lij are the onsgger coefficients. The Onsager equation for the flux of heat Jq is -Jq = 2 Iqi[aEi/ax]T + qu[alnT/dx]. (2-3) The matrix of Onsager coefficients is symmetric. lij a ljiv i.j = 1.2. (2-4) i = 1.2. ‘ (2-5) The electrochemical potential of ion i is given by n. = u“: + 1n[y.c.] + 2.1% (2-6) where yi is the molar activity coefficient. zi is the charge number of ion i. F is Faradays constant. and 6 is the electrostatic potential. Dudley and Steele [1980] and Horne [1981] write the gradient of ii in terms of the independent °i and 6 as aEi/ax = E pik(3ck/ax) + ziF(36/dx) — si(aT/ax) (2-7) 8 where Si is the partial molar entropy of ion i with [aEi/axjr - [aEi/ax] + si[aT/ax]. (2-8) uij . [aEi/acj] . (RT/cj)[bij + (alnyi/alncj)], (2-9) where 5ij is the Kronecker delta. With the substitution of Eq. (2-7) into Eq. (2-2). the diffusion fluxes become for each ion -Ji . 2 Dij[acj/ax] + (Ati/in)[ad/6x] + liq[alnT/ax]. (2-10) with Dij a 2 1m“. ti - [ziFZIMJ 2 zklik. A = F222 zizjlij' Eti = 1. (2-11) The requirement of reciprocity imposes upon the diffusion coefficients the relation D12911 ‘ D11ll12 = D211122 “ D22F21- (2'12) Eq. (2-9) reduces to the Nernst-Planck equation under the necessary and sufficient conditions that (Horne [1981]) i) V030, 192. (2-13) H. N iii) ln(yi) . O, The Nernst-Planck equations are valid only for ideal solutions. in which moreover there is no diffusional coupling. The condition of local charge neutrality implies that the current of one species must everywhere be counterbalanced by that of the oppositely charged species. This coupled motion can be described by only 532 independent velocities (electrolyte and solvent). and accordingly only 233 independent diffusion flux. For nonisothermal diffusion we have as well only one independent thermal diffusion coefficient. Haase [1969] defines the diffusion fluxes under these conditions as ’J ’ D[(ac/ax) - oc(aT/ax)] (2-14) where c is the concentration of the neutral electrolyte. o is the Soret coefficient. and D is the diffusion coefficient (also known as the "ambipolar" or "chemical" or "mutual" diffusion coefficient). The diffusion coefficient defined by Eq. (2-14) can be related to the Di” J of Eq. (2-9) by D ' 2(D11D22-D12D21)/(Dll + D22 ' (21/22)Dlz ’ (22/21)021) (2’15) Dudley and Steele [1980] evaluate this diffusion coefficient in terms of 10 what they call a) thermodynamic factor. Wagner [1975] derives an expression for this same diffusion coefficient for the case of a metal excess 8 diffusing in compounds of stoichiometric composition A1+va in terms of the Onsager coefficients and chemical potential gradients. Horne [1981] directly relates the diffusion coefficients of Eq. (2-9) to the mutual diffusion coefficient D for a symmetrical electrolyte. For a nonisothermal system the Soret coefficient a in Eq. (2-14) can be related to the individual ionic heats of transport by c - -O‘/Tc(ap/ac) (2-16) where Q. is the heat of transport of the electrolyte given by Q. = leI + V205 and v1 and V; are the stoichiometric coefficients of species 1 and 2 respectively. The individual ionic heats of transport are defined by ‘ # 1iq ‘ 2 1ij°j~ (2-17) With Eq. (2-17) the heat flux can be rewritten in terms of the heats of transport as - H P _ Jq a -2 Ii Oi + x.(aT/ax) (2 18) where KS a qu/T — 2 (lqu;)/T. (2-19) 11 C. Continuity Equations In addition to the flux equations of nonequilibrium thermodynamics given above. the system will also obey the relevant conservation of mass equations given by (aci/at) = - (OJi/ax) + 2 E virbr (2-20) where t is time. vir is the stoichiometric coefficient of species i in chemical reaction r. and hr is the rate of chemical reaction r. An additional transport equation comes from the conservation of energy which for a nonisothermal system is written as Haase [1969] OC§(aT/6t) 3 “(qu/ax) - 1(36/ax) "2 2 Virbrai - 2 Ji(aEi/ax) (2-21) where p is the system density. C p is the specific heat. I is the total current. and Hi is the partial molar enthalpy of component i. An excellent approximation to this equation for problems in which there are neither reactions nor current.flow (Horne and Anderson [1968]) is pC§(aT/at) = (a/ax)[x.(aT/ax)]. (2-22) D. Electrostatic Equations 12 A fundamental law of electromagnetism (Bleaney [1957]) which relates the electrostatic potential to the concentrations of charged species in the system is Poisson's equation and is given by e(aE/ax) = F 2 ciz (2-23) is where a. the permittivity of the medium. is considered constant. and the electric field E is defined by E = -(36/ax). (2-24) Taking the time derivative of the left hand side of Poisson's equation. and substituting Eq. (2-20) for the right hand side. we find e(a/at)(aE/ax) = -F§ zi(aJi/ax), (2-25) where we have also used 2 virbrzi = 0. (2-26) If the electric field and charge density as well as their space and time derivatives are continuous functions. then it is permissible to interchange the order of differentiation in Eq. (2-25). giving (a/ax)[s§ ziJi + s(8E/6t)] = 0. (2-27) 13 Ampere's law requires that aI/ax = 0. (2-28) and therefore it would be incorrect to define the total current as being equal to the summation term in Eq. (2-27). which accounts only for the motion of the electric charges. because this would conflict with the continuity equations whose validity is confirmed by all experiments. Naxwell [1865] was the first to reconcile this difficulty and correctly defined the current as I = F 2 11-71 + MOE/at). (2-29) where the second term. the Maxwell displacement current, arises because of the time dependence of the electric field. E. Combined System of Equations The displacement current equation. Eq. (2-29). along with the conservation of mass and energy equations coupled with the nonequilibrium thermodynamic equations provide the basis for a complete macroscOpic description of the system in terms of material properties and the dependent V'ti'bIOS 01. c2. and E. The full set of equations is obtained by substitution of Eq. (2-10) into Eqs. (2-20) and (2-29): (aci/at) a (a/ax)[ 2 Dij(acj/ax) + liq(alnT/3x) 14 - (Ati/in)E ] + E virbr. i=l.2. (2-30) e(aE/aT) = I + F 2 2 ziDij(8cj/8x) + F 2 ziliq(dlnT/ax) - IE. i.j=1.2 (2-31) where the solvent reference velocity has been neglected. Leckey and Horne [1981] distinguished electrically neutral behavior such as ambipolar or mutual diffusion from behavior due to non-electroneutrality and introduced the sum 8. and difference or charge density A composition variables. U.) I - (1/2)[(c2/z1) - (cl/22)]. (2-32) > I - (1’2)[(c2/21) + (cl/12)]. . (2-33) The transport equations become aS/at = (a/ax)[ Tss(33/ax) + YSA(3A/ax) + qu(3T/6x) - (8Y33/2z122F)E ] + (1/22122) 2 (vzrzz - v1,11)b,. (2-34) aA/at = a/ax[ YAs(aS/ax) + YAA(8A/8x) + YAq(8T/8x) q - (aYAE/ZzlzzF)E J, (2-35) 15 3(aE/at) = I — (IIZzlzzF)[ YAs(aS/6x) + YAA(aA/8x) + YAq(aT/ax) - (eYAE/2zlzzF)E J, (2-36) where thfl Yij (i.j = S.A.E.q) are defined in Appendix A. These transformed equations highlight the very different physical behaviors associated with diffusion of neutral matter on the one hand. and diffusion of charged matter on the other. Eqs. (2-34) through (2-36) along with the conservation of energy given by Eq. (2-21) form the complete set of partial differential equations whose solution we seek. The solution will depend upon the system's material parameters and the boundary conditions placed upon it. F. Solving the Equations The nonlinearity of Eq. (2-21) and Eqs. (2-34) through (2-36) makes it highly unlikely that one can obtain analytical solutions to them in .closed form. The Yij are in general concentration and temperature dependent and so is the conductivity. For systems in which there are reactions some sort of model dependent behavior for the reaction rate must be postulated. Additionally much information is required in order to define the material parameters necessary for expressing a solution. For solid state systems such information is 808100. In addition to attempting analytical solutions to the transport 16 equations. it is necessary to Obtain numerical or computer solutions as well in lieu of analytical solutions. Cohen and Cooley [1965] were the first to simulate ion transport on a large scale. Several others (Kahn and Naycock [1966]. Feldberg [1970], Sandifer and Buck [1975]) have numerically solved the simplified set of transport equations (Nernst-Planck equations) under a variety of conditions. Brumlevé and Buck [1978] have developed the most general program of this type. The prOgram used in the research which led to this dissertation was developed by Leckey [1981]. and extends the generality of the Brumlevé and Buck program by using the full set of transport equations (Eqs. (2-34) to (2-36)) valid at any concentratiOn. for isothermal systems of uni-univalent electrolytes. Here. the transformation of the full set of transport equations to finite difference form is extended to include both electrolytes of any valence and nonisothermal systems. The use of the full set of equations allows one the option of incorporating cross effects and activity coefficient effects. both of which are ignored in the Nernst-Planck formulation. A listing of the complete set of dimensionless variables chosen for this system is included in Appendix C. A listing of the program "IONFLO" is shown in Appendix F. A more complete discussion of the details surrounding the transformation of the transport equations to finite difference form can be found in Chapter 6. The principal value in carrying out numerical solutions to problems like this one is that one can (1) model the behavior of various material parameters based on available data. and (2) test models by comparison 17 with experiment. Numerical solutions also provide a check on the accuracy of analytical solutions. they allow verification of the assumptions involved in obtaining them. and they provide clues to improved analytical solutions. CHAPTER 3 BOUNDARY CONDITIONS A. Introduction For the one dimensional system described in Chapter 2 there are two boundaries. At each boundary it is necessary to make a statement about the concentrations of the ions. and for nonisothermal systems it is necessary to make a statement about the temperature as well. In order to make contact. with experimental results boundary conditions must be devised which enable one to simulate as closely as possible actual experimental conditions. For applications considered here. battery-type devices, boundary conditions of the type Ii = kifcil ‘ kibcio (3-1) previously used by Brumlevé and Buck [1978] are employed. The flux of an ion i at an interface is described by a forward and reverse rate constant kif and kib and ion concentrations in each phase °il and °i0° The direction of flux and interfacial concentrations are shown in Figure (3-1). The assumptions required for the validity of boundary conditions of 18 19 Figure (3-1). Interfacial region at a phase boundary. ELECTRODE 20 INTERFACE kif + kib J1 +-—————5——————-+ = k C if 11 'k Figure (3—1) ELECTROLYTE 21 this type have been discussed previously by Horne. Leckey. and Perram [1983]. For closed systems application of the above boundary conditions is a simple matter. in that all interfacial rate constants are set equal to zero. and no matter traverses the boundaries of the system. For open systems where electrode reactions may be occurring it becomes necessary to assign to these constants "reasonable" values that reflect the behavior of 'the particular electrode system. This requires additional experimental data in the form of measurements of the a.c. and d.c. prOperties of the electrode-electrolyte system of interest. Braunshtein. Tannhauser. and Riess [1981] have measured these properties for platinum electrode-doped metal oxide systems. In the event that the electrode system behaves in a completely reversible manner. Ji I 0 and the rate constants in Eq. (3-1) approach infinity while °i0 9 °il (Leckey [1981]). The application of these boundary conditions becomes more difficult for real systems (k's # 0.”) and requires that we look at the physical interpretation of the terms in Eq. (3-1). B. Physical Interpretation of the Boundary Conditions The three assumptions which lead to the kinetic flux expressions given by Eq. (3-1) (Horne. Leckey. and Perram [1981]) are: i. The cross terms within the interfacial region are negligible. ii- The quantity “3 + ziFd is a linear function of x through the interface. iii. The flux Ji and diffusion coefficient Di are constants through the interface. These assumptions result in the following expressions for the rate 22 constants. kif = (Di/8) [ aiexp(-ui)/(l-exp(-ai)) ]. (3-2) kib = (Di/8) [ ai/(l-exp(-ai)) ], (3—3) where “i = AI»? + ziFdllRT. (3-4) and 8 is the interfacial thickness. The standard state chemical potential is included here because it is a property of the phase in question and in general we are speaking of the exchange of matter between two different phases. For phase boundaries at which [ziFAOI )) [Aug]. a large positive value of “i corresponds to the situation in which a positive ion is moving to a region of more positive electrostatic potential. or a negative ion is moving to a region of more negative electrostatic potential. In either event the magnitude of the rate constant for movement in this direction (kif) is expected to be small while the rate constant for movement in the reverse direction (kib) should be large (see Figure (3-2)). Conversely. when “i is a large negative number exactly the Opposite should be the case. namely movement of positive ions to a region of more negative potential. or negative ions to a region of more positive potential. In either case we e31,99t th‘t kif becomes large. while kib becomes small (see Figure (3‘2)). In cases where [Aug] >> [ziFA6|. Eqs. (3-2) and (3-3) predict movement of ions to the phase in which their standard state chemical 23 Figure (3-2). Dependence of the kinetic rate constants on a1. 24 ANIMV unawam . E>e<..__~ + .33 cod cod oofil _ + u + + 1 A I F 8... L1 11 8., 1... 11 SN LI IT 8.,” L1 11 8... a; f A 1 .1 A A 1 .i + .. 8S (lg/9) x lUDlSUOQ 3303 25 potential is lower. (Aug) -) -.. implies ai -> .. and hi. -) o. Perram [1980] approximates the value of A“? from the Born equation An‘i’ - [(ziF)2/2si] 131/er) — (mm (3-5) where 'i is the radius of ion i. and er and ‘1 are the permittivities of the right and left hand phases respectively. In accordance with the argument presented above. when 81 < gr, A“? ( 0, and transfer of material to the phase with the higher dielectric constant is favored. The application of the kinetic boundary conditions given by Eq. (3-1) should simulate as closely as possible the experimental conditions. Under certain conditions it is possible to make quantitative statements concerning the values these constants should assume. For instance when an electrode is blocking to a particular ion. that ion cannot penetrate thefiinterface from either direction, and one can assign values of zero to the rate constants in Eq. (3-1). From Eqs. (3-1). (3-2). and (3-3) we see that this can be true if and only At the other extreme there is the reversible electrode at which Ji 8 0. implying kif°il 8 kibciO. In cases where “i is zero we see from Eqs. (3‘1) and (3-2) that kif a kib' and therefore cil B ciO. In cases where neither reversible nor blocking conditions exist at the boundaries. the matter of assigning reasonable values to the rate constants is complicated by the fact that it is necessary to know oi in 26 Eqs. (3-2) and (3-3) in order to obtain values for the rate constants. The change in the value of “i from its equilibrium value of zero to its nonequilibrium value is due entirely to the change in the value of the contact potential (A4 in Eq. (3-4)). since the value of Aug remains constant regardless of current flow. The change in the value of a is defined as the overpotential. and is not experimentally measurable. For two phase junctions in series (A01; 4' Ac§)/ziF = 11L + “R. (3'6) where the overpotentials ("L and DR) are “L a (A‘ _ A‘rov)Ls (3-7) and "R = (A6 - A‘rev)R' (3-8) The total overpotential loss (“L + UR) for a given electrode electrolyte system can be estimated from the analysis of the impedance frequency response of the system. This analysis is often aided by analogy with equivalent electrical circuits consisting of resistors and capacitors each associated with a single process (Sandifer and Buck [1975]. Warburg [1899. 1901]). In Figure (3-3a). the bulk resistance is represented by Rb. while the resistance at the interface is represented by Rig. The charging of the double layers through the surface resistances is associated with surface capacitance Cif° Figures (3-3b). and (3-3c) show the equivalent circuit representations of ideally 27 polarized and depolarized electrode-electrolyte interfaces respectively. Impedance frequency response data (Braunshtein. Tannhauser and Riess [1981]) for a given system can then be used to assign values to the resistances and capacitances found in an equivalent circuit representation. which leads to estimates of the magnitude of the total overpotential given by Eq. (3-6). C. Applications of the Boundary Conditions Various workers have used these boundary conditions to simulate a variety of conditions. For liquid-liquid junction cells with blocking electrodes Leckey [1981] set the rate constants equal to zero at each end. .He also used these boundary conditions to study charge injection into a one dimensional cell at a blocked interface containing a reversible electrode at the other boundary: For the reversible electrode simulation the rate constants were set. to machine infinity. Chang and Jaffe [1952] studied the effect of changing the values of the rate constants on the conductance. Brumlevé and Buck [1978] simulated a variety of effects ranging from charge steps at blocked interfaces (k's = 0) to ion exchange membranes (k's finite). Their assignments for the rate constants were made from analyses of the impedance frequency responses of the systems as described above. For the one dimensional solid state systems discussed in Chapters 4 and 5. it is convenient to break up the electrolyte into three regions (Figure (3-4)). In the vicinity of each electrode (regions I and II) there are reactions occurring, while in the bulk between regions I and 28 Figure (3-3). Circuit representations of solid-solid interfaces. (a) Partially blocking electrode. (b) ideally polarized electrode. (c) ideally depolarized electrode. 29 b Bulk Bulk Interface l I Interface Rif Bulk Figure (3-3) Interface 30 Figure (3-4). Subdivision of the electrolyte into three pseudophases with pseudophase boundaries shown at x = t L/2. 31 L ,R kif kif ___.————-+ ————-——+ L R kib kib +——-—— e-—————— ANODE (I) ‘ (C) II CATHODE ' Bulk Electrolyte _ ‘ -L/2 . +L/2 ~ Figure (3-4) 32 II there are no reactions and only transport is occurring. The lines that separate the reaction regions from the bulk denote pseudophase boundaries. The term pseudophase is used here because there is no discontinuity in system prOperties at these boundaries. only a difference in the conservation of mass equations on either side with aci/at = -(BJi/ax) (3-9) in the bulk and aci/at = -(81i/ax) + 2 Virbr (3-10) in regions I and II. Therefore at the pseudophase boundaries An? = 0 and Ad ' 0. and by Eq. (3-4) “i s 0. From Eqs. (3-1). (3-2). (3-3). and (3-10) the flux of ion i at each boundary becomes Ii 3 (Di/5)(Cil ’ 01°). (3’11) 1‘if = l‘ib = 131/5. (3-12) Eq. (3-11) applies whether charge transfer between the electrode and electrolyte is reversible or irreversible. When the reactions are at equilibrium. b, = 0 in Eq. (3-10). and Ji(Bulk) = Ji(I) = Ii(II) = 0 with cil = °i0 by Eq. (3-11). In general the manner in which the component fluxes in the reaction regions change from their bulk solution values to their values at the 33 electrode electrolyte interface is determined by the behavior of the electrode (overpotentials) and by the magnitude of the reaction rates in that region. At the steady state Eq. (3-10) gives J1(1 or II)) = I 2 (virbr)dx. (3—13) where the integral is taken from the electrode to the bulk solution. Knowledge of the dependence of the reaction rates on distance from the electrode would allow one to obtain an analytical solution to Eq. (3-13). where the behavior of the electrode towards species i determines the integration constant. High values for the rate constants used in Eqs. (3-11) and (3-12) imply electrodes behaving nearly reversibly with very little change in ionic concentrations between bulk and reaction phases. and hence very little double layer charging. From the equivalent circuit point of view the value of the interfacial rate constant is a measure of the resistance of the double layer. High rate constants mean low interfacial resistances and little double layer charging. Low values of the rate constants imply that significant concentration differences will arise across the pseudOphase boundaries and double layer charging will occur. _These rate constants correspond to high values of interfacial resistances. In Chapters 4 and 6 the effect of charging the double layer is examined through inclusion of the charge density terms in the boundary condition formulation. CHAPTER 4 STEADY STATE TRANSPORT IN OPEN SYSTEMS A. Introduction In this chapter the equations of nonequilibrium thermodynamics are Used to investigate the steady state behavior of energy conversion devices (batteries and fuel cells). Of special interest are solid state electrolytes into which a mobile Species can enter from the static lattice of the electrodes. Intercalation electrodes (Whittingham [1979]) are useful for this purpose since they provide near reversibility of reaction. The equations are. however, quite general and can be applied with appropriate modification to other electrode and electrolyte systems. B. Description of the system The fuel cell arrangement (Figure 4-1) (Kiukkola and Wagner [1957]) consists of a solid electrolyte (C) sandwiched between two electrodes (B.B') not necessarily at the same temperature. Each electrode is attached to a wire (A.A') which leads to a terminal at some mean temperature TA- The measurable EMF is A6 = ‘VIII - 61. The model electrolyte (Figure 3-4) extends from x = -L/2 to x = +L/2 and the voltage drOp across it is 6v - ‘IV (Figure 4-1). In general there are 34 35 Figure (4-1). Nonisothermal galvanic cell. Wire (A) 36 TB TB, Anode ' Electrolyte Cathode (B) (C) (B') ¢III ¢IV ¢V ¢VI Figure (4-1) BI ¢VII Wire (A') VIII 37 reactions occurring near the electrode electrolyte interfaces which provide the driving forces (gradients of chemical potential and temperature) that cause the flux of material through the electrolye to occur. We assume in general a one dimensional problem with the goal of determining the voltage loss across the electrolyte (JV-61v) using the equations of nonequilibrium thermodynamics. This voltage loss can then be subtracted from the theoretical or Nernst potential. and an estimate of the thermodynamic efficiency of such a device can be obtained. This estimate neglects overpotentials. which are changes in heterogeneous or contact potentials (‘VII - ‘VI' “VT - 6v. ‘IV - ‘III' and 6111 - 611) from their equilibrium values. Contact potentials arise as a consequence of charge transfer from one phase to another. Estimates of these contributions to the total voltage loss across any Operating device can be obtained from the measurement of both the a.c. and d.c. prOperties of the system of interest. Such measurements lead to improved estimates of the true operating efficiency of such devices. C. Steady State Equations At the steady state the time derivatives of all intensive variables vanish and the fluxes defined by Eqs. (2-10) reduce to -J"{ - E Dij(dcjldx) - (Ati/ziF)E + liq(dlnT/dx). (4-1) where I: represents the constant steady state value of the flux of ion i. The sum and difference pseudo fluxes are defined in the same way as the sum and difference concentrations, 38 IX - (1/2) [(13/21) + (If/22)]. (4-2) I; = (1/2) [(13/21) - (If/22)]. (4—3) From Eqs. (2-34) through (2-36) we see that in terms of dependent variables and material parameters. -J§ . Yss(dS/dx) + Y5A(dA/dx) - YSEK'IB + qu(dT/dx). (4-4) -JX - YAS(dS/dx) + YAA(dA/dx) - YAEx’lE + YAq(dT/dx). (4-5) where K = (23122F)/e. (4—6) The pseudo fluxes IE and I: are not true Onsager fluxes since there is no reciprocal relationship involving only the Yij Of Eqs. (4-4) and (4-5). Simplification of these equations results from simultaneously solving them for (dS/dx) and (dA/dx): dS/dx = as + BSE + 13(dT/dx). (4-7) dA/dx = “A + 3A5 + 7A(dT/dx). (4-8) where the coefficients a. B. and y are defined in terms of the Y-. in 1J 39 Appendix A. In general these coefficients incorporate all effects due to nonideality of the system and are themselves dependent upon the concentrations of all species present as well as the temperature. It is this dependence that gives rise to the nonlinearity of Eqs. (4-7) and (4-8). The a's are related to the fluxes at the boundaries and would be zero for the case of the system sandwiched between blocking electrodes. For ideal systems for which the Nernst Planck equations are valid the c's behave as constants. while the 3's are linear functions of the ionic concentrations. The 1's incorporate nonisothermal effects such as the heats of transport of the ionic species. For systems in which there are small temperature gradients these coefficients are nearly constant. with values dependent upon the mean temperature of the system (Cowen [1979]). It is possible to eliminate A as a variable by combining Eq. (4-8) with Poisson's equation to give (A’s/4x2) - (KBA) - (KaA) + K1A(dT/dx). (4—9) Eqs. (4-7) and (4-9) are the starting set of equations for the macroscopic investigation of the steady state problem. The only assumptions made up to this point are those implicit in the Onsager formulation of the transport equations of nonequilibrium thermodynamics (Haase [1969]. Horne [1966]. and Onsager [1931]). All effects due to nonideality are retained. D. Solution to the isothermal problem 40 For the case of no temperature gradient. the nonisothermal terms in Eqs. (4-7) and (4-9) vanish. giving 2 2 _ (d E/dx ) - (KBA)E - KOA (4 10) dS/dx 8 as + SSE. (4-11) For ideal solutions Eq. (4-10) reduces to the Poisson-Boltzmann equation with KBA the reciprocal of the Debye length squared. KBA - (2F2I°)/SRT. (4—12) and Ic is the ionic strength of the solution. The boundary conditions for the Open system are I; de - Constant. (4-13) A = 9L. 1 = -L/2. (4—14) A ‘ PR) x '3 +L/29 (4-15) where 9L and 9R are the charge densities of the system at the left and right hand walls respectively. Their values will depend upon the speed of charge transfer across the interface separating the Open system from the adjoining phase. The first boundary condition is simply a statement that the total mass of the system is fixed at the steady state. 41 Because of the nonlinearity of Eqs. (4-10) and (4-11). it is generally not possible to obtain solutions to them in closed form. For this reason a pertubation approach is used in which system parameters and variables are characterized according to their departures from their mean values. Before introducing this technique it is useful to transform the starting equations to dimensionless form in order to facilitate comparison with numerical results. The pertubation scheme is simplified simultaneously. The following dimensionless variables are introduced s - (s - s.’/S.' aj - (AA/Sn)cj. j . s.A. ‘- - _ 2 A - A/sn. pj - [( 221z2F(AA) )/.o]aj. j - S.A. x - x/AA. 'E . -[eo/(22122F8n(AA))]E. (4-16) where S. is the mean value of the sum concentration. so is the permittivity of a vacuum. and the barred quantities represent dimensionless variables. The value of AA in the above equations is chosen such that the reduced Debye length is fixed at a value of 5.0. This number has been found to give a convenient grid spacing length for numerical simulations (Brumlevc and Buck [1978]). The dimensionless variables used here are the same as those used in subsequent numerical simulations. In terms of reduced variables the starting equations become (dz-EM?) - (Pi - —;A/EP_S. (4-17) 42 (dS/d'x') - Es? + ESE. (4-18) where EPS is the dielectric constant of the system (78.5 for water) and X is the reciprocal of the reduced Debye length of the system. given by I - JEAI'EPS (4-19) In reduced variables the boundary conditions become I; 3‘13 " 0: (4-20) I - 91/3... I =- —L/2AA. (4-21) I - .313... 2 = +L/2AA. (4-22) The pertubation scheme is defined for the dependent variables as follows 8 = so + 881 + 682 + °'°, (4-23) E-Eo+ael+aaz+ (4-24) where 5 is a bookkeeping device that allows one to characterize the departure of Eqs. (4-17) and (4-18) from linearity. The barred notation has been omitted from Eqs. (4-23) and (4-24) for simplicity. The coefficients (L) are expanded in a Taylor's series about their mean 43 values according to L(S.A) = Lo + 8[LS(S) + LA(A)] + °-°. (4—25) where L0 represents the value of L at the mean values of S and A. and Lj (j = S. A) is the derivative of L with respect to j evaluated at the mean values of S and A. Substitution of Eqs. (4-23) through (4-25) into Eqs. (4-17) and (4-18) and subsequent ordering of the terms according to their power in 6 results in sets of coupled linear ordinary differential equations solvable by standard techniques. Although in principle this set of equations is solvable. there is no guarantee of success in this technique. as the magnitude of higher order terms may increase leading to a divergent solution. Collection of terms of order zero in 8 results in (donldxz) - AfiEo = -aA0/EPS. (4-26) (dSoldx) 8 “$0 + BSOEO' (4-27) where Lij is the jth order approximation to the coefficient Li (1 = S,A, j ' 0.1.2.°°'). This system of equations would be the true description of the system if it were ideal. To obtain E0. we first solve Eq. (4-26). To obtain So. we use the 44 solution to Eq. (4-26) in Eq. (4-27). The zeroth order contribution to the charge density (A0) is obtained by simple differentiation of E0 in accordance with Poisson's equation. The general solution shown below results in the generation of three arbitrary constants whose values are fixed by application of the boundary conditions. 'e require that the zeroth order solution satisfy totally the charge density requirement at each wall while higher order solutions in charge density will be zero at the walls. The boundary conditions then become I: Sidx . o. i - o.1.2,--~, (4—23) A0 3 pL/Sn. x a -L/2AA. A0 '- pR/Sn. x = +L/2AA. A1 = 0. i - 1.2.-~-, E - i(L/2AA). (4-29) The zeroth order solution is 30 - -(qu/BA0) - (p/A(EPS)) [cosh(Ax)/sinh(A§)]. s0 . (szo/BAO)(sinh(Ax)/sinh(A§)) + [(asoaAo - aAoBso)/9Ao]" D0 - p [sinh(Ax)/sinh(A§)]. (4-30) where the dimensionless cell length defined by 45 E ' (L/2AA) ~ (4-31) and p the charge density p = “(PL/3..) - (pk/Sn) (4—32) have been introduced. It is apparent from the form of the solution that the redefinition of the concentration variables has served an exceedingly useful purpose in that each of the first two expressions is separated into two parts. one due to charge density and one due to bulk or electrically neutral behavior. At 298 K. for an aqueous solution of ionic strength 0.01M. and a cell length of 1 cm. the value of 5 is on the order of 108. For solid electrolytes (where dielectric constants are much lower) at temperatures up to 1000 K and at cell lengths of up to 1 cm. the value of t is on the order of 108 or higher. This means that contributions to E0 and SO due to nonelectroneutrality are going to be negligible everywhere except within approximately 5 Debye lengths of the walls. In the bulk solution a linear drop in concentration with a constant value for the electric field is found as expected. Collection of terms of order one in 8 gives (d2E1/dx2) - (A30)E1 - -(1/EPS)[BASSO + BAAAO]EO. (4-33) 46 (dslldX) ' 33031 I [53380 + BSAAO]EO' (4‘34) where the Lij' (i.j - S.A) are the derivatives of Li with respect to j evaluated at the mean concentrations of the system. Inclusion of terms involving derivatives of the u's is deferred until higher orders because of their expected smallness. The procedure for solving the first order set of equations is similar to that for the zeroth order. first the expressions for So and B0 are substituted into Eq. (4-33). and a solution for 81 is obtained. Then. E1. along with S0 and E0 are substituted into Eq. (4-34) to obtain 81' As before A1 results from simple differentiation of E1 in accordance with Poisson's equation. The complete first order result is shown in Appendix B. Higher order solutions are obtained in the same way. The inclusion of more and more terms in the starting equations results in more and more complicated algebraic problems as one progresses to higher orders. This is already apparent in the first order case. The way the problem was defined and the nature of the pertubation scheme result in solutions that alternate in parity. For the zeroth order case. E0 is an even function while So and A0 are odd; at first order. 81 is odd while 81 and A1 are even. This requires that one obtain a second order solution for the electric field in order to determine the first correction to the voltage drop Ad across the electrolyte since. for E1 odd. A61 - - I: Eldx = o. ‘ (4-35) 47 The results are shown in Appendix B. E. Determination of Efficiency The current delivered to the load is I 8 F 2 2111 I 11 + I2. where 11 and I2 are the individual ionic currents. The efficiency a defined here is the product of the voltage efficiency (A¢)/(A¢ ) and the rev faradaio efficiency (I/In). where In is the maximum current that would be produced by the cell if it could operate ideally. The power P delivered to the load by an ideally Operating fuel cell (Bockris [1969]) is P = 1A6 = I[A¢rev - IL/Am], (4-36) where Am is the mean specific conductance of the electrolyte. IIn is easily determined by setting the derivative of the power with respect to the current equal to zero. giving Im = AmAerev/ZL. (4—37) Measurement of the a.c. end d.c. prOperties of the electrode- electrolyte system of interest enables one to determine the value of Am and to calculate faradaio efficiencies as a function of total current produced by the cell. The electrochemical contribution 6v - ‘IV to the total measured 48 voltage across the cell A6 is obtained from integration of Eq. (4-30). In the isothermal case. TA = Th 8 Th: so that 611 = 61. and ‘VII 3 6v111 (see Figure (4-1)). The measured EMF consists of the sum of the heterOgeneous contributions A‘het and the homogeneous contributions A‘hoxav A‘ ‘ A‘hom + A‘het = AVIII ‘ 61. (4‘33) where Aehom = 6V - ‘IV' (4-39) A‘het = UVIII ‘ ‘VI’ + (‘VI ‘ Av) + (‘Iv ' “111’ + (5111 - 51), (4—40) The assumption of electrochemical equilibrium across the wire/electrode interfaces allows us to write fie(A') = Ee(B'). fie(B) = He(A). i,(B) = He(x=-L/2). fie(3') = He(x=+L/2). (4—41) By substituting the full expression given by Eq. (2-6) for the electrochemical potentials above and solving for A‘het' we obtain FAdhct = (RT)ln[a,(-L/2)/ae(+L/2)]. (4-42) where ac. the activity of the electron, is different at either end of 49 the electrolyte. Because both the standard state chemical potential and the activity of the electron are equal in phases A and A' since the measurement of A6 must take place between wires of the same phase. they do not appear in the final expression for the heterogeneous potential. The heterogeneous contribution is then added to the homogeneous contribution. FAd = F(6v-s1v) + (RT)ln[ae(-L/2)/ae(+L/2)], (4-43) which enables one to determine the voltage efficiency of the Operating device. From a knowledge of the total efficiency of the Operating device as a function of current given by the macrosCOpic description of the system, it becomes possible to predict under what external loads the cell operates most efficiently and how the material parameters affect cell performance. F. Results and Discussion 1. Introduction The results and ensuing discussion presented in this section are based upon the prOperties of the mixed (ionic-electronic) conductor solid electrolyte doped Cerium (1V) oxide. Choudhury and Patterson [1971]. Tannhauser. [1978]. and Riess [1981] have all derived analytical expressions relating the steady state voltage output of high temperature (T 3 1100 K) fuel cells employing this electrolyte to the current produced by the cell. Such relationships are termed "characteristics." 50 and depend upon the material parameters of the electrolyte. Their calculations all share the common shortcomings of using the Nernst-Planck transport equations rather than the full set of Onsager flux equations. and at some point in their calculations all assume electroneutrality throughout the electrolyte. In addition they neglect the spatial dependence of the various system parameters (D. I. and K) and in effect consider only the zeroth order solution for an ideally behaving system. These assumptions could. collectively. lead to a serious overestimation of the efficiency and hence of the economic performance of such systems. From the data acquired by the above workers and the more accurate expressions for the potential derived here it is possible to assess the importance of these effects and to arrive at a more realistic theoretical prediction of the economic feasibility of such devices. 2. The Mixed Conductor System The electrolyte is a membrane of doped CeOz, typically at high temperature. extending from x=-L/2 to x=+L/2. The arrangement and phase boundary sections are shown in Figures (3-1) and (3-4). The doping is done with a lower-valent metal. Ca2+, creating vacancies V, in the oxide sublattice with an effective double positive charge above what exists in the undoped material. These charged vacancies conduct the ionic portion of the electrical current. The overall reaction which provides the driving force for the cell is given by the sum of the reactions occurring at the phase boundaries, 51 “.203 + GFCO - “.3040 (4‘44) The oxygen partial pressures at the phase boundaries are fixed by the equilibria in the two phase electrode compartments. Anode: 6Fe0 + 05(g) = 2Fe304. (4-45) Cathode: 6Fe203 - 4Fe304 + 0§(g). (4-46) Substitution of Eqs. (4-45) and (4-46) into Eq. (4-44) produces 02(s.B) - 02(s.B') (4—47) with a AG. equal to that of the reaction given in Eq. (4-44). This value of AG. is the amount of reversible work required to displace one mole of 02(g) from its partial pressure at B' to its partial pressure at B. The reversible potential is therefore A"rev ' (RI/“F, ln [PO2(g.B')/P02(g.8)]' (I-48) where n is the number of equivalents of electrons required to displace one mole of 02. and P02 is the partial pressure of oxygen. Due to the kinetic problems at the phase boundaries and the ohmic losses that occur in the electrolyte the useful voltage obtainable from such a cell will always be less than Aerev. For the reactions given by Eqs. (4-45) and (4-46) the equilibrium 52 conditions at either end of the electrolyte are (Schmalzried [1974]) phase I: 02— - 2e- + (1/2)02(g) + V3. (4-49) phase II: v, + (1/2)02(a) + 2a” a 02'. (4-50) Since the partial pressures of oxygen are fixed in the two phase electrode compartments. the electrochemical potentials of the electrons. vacancies. and oxide anions will in general be different at each end of the electrolyte. and therefore a flux of each species occurs. The entropy production (9) for this system is (Howard and Lidiard [1963]) e = Jv'xv' + J. x.- + Joz-x02—. (4—51) where xi - - (ail/ax). , (4-52) The conservation of lattice sites requires that 102- = -Iv'. since every jump of an oxide anion from one lattice site to another is accompanied by a vacancy jump in.the opposite direction. Therefore 0 becomes -Jv, [ (apv./ax)-(au02-/ax) ] - Je- (duo-lax). (4-53) and Onsager thermodynamics yields -11 - 111 [ (apl/ax) — (auoz-Iax) - 202-F(66/ax) ] 53 + 112 (apzlax). (4-54) -12 = 121 [ (Bel/6x) - (BuOZ-Iax) — zoz-F(ad/ax) ] + 122 (Buzlax), (4‘55) where the subscripts 1 and 2 denote vacancies and electrons respectively. and Eq. (2-7) has been used. Because the number of vacancies is far less than the number of oxide anions. the relative change in the oxide anion concentration will be far less than that of the vacancies. We therefore neglect the change in the oxide anion chemical potential. Eqs. (4-54) and (4-55) simplify to the prototype flux equations for isothermal tranSport given by Eqs. (2-2) -Ji = 2 Iij[aHj/axj. i.j = 1.2 (4-56) where 21 8 202- is the effective charge on the vacancies. 3. Calculation of Efficiency From the definition of efficiency and Eqs. (4-37) and (4-43) we write n = [21L/Am(A6rev)2] [ (av-61v) + (RT/F) .ln[cz(-L/2)/c2(+L/2)] ] (4-57) 54 where the activities have been replaced with concentrations. These concentrations subsequently serve as boundary conditions through the equilibrium constants for the reactions given by Eqs. (4-49) and (4-50). The value of ‘VP‘IV_in Eq. (4-57) is determined by integration of Eq. (4-30) and allows us to determine the effect of selectively changing the value of p on the efficiency through Eq. (4-57). All effects due to the nonideality of the electrolyte are retained in the evaluation of ev-JIV. The more mobile electrons enter and leave the electrolyte faster and therefore create a region of negative charge at the anodic end of the electrolyte and a positive charge at the cathodic end. Integration of Eq. (4-30) reveals a linear relationship between the yoltage loss and this charge buildup. (p). (av-in) - . + bp. _ (4-58) with the values of a and b displayed as a function of current in Table (4-1). These data highlight the importance of the inclusion of space charge upon voltage loss estimates. Only when one approaches complete charge separation (p 9 1) does this effect become important. but it is precisely at the charged surface of the electrode where one might expect complete separation of charge to occur. From Eq. (4-57) a linear relationship between the efficiency and the charge density results with n = c + dp. (4-59) 55 The values of the coefficients c and d for Eq. (4-59) are tabulated alongside those of Eq. (4-58) in Table (4-1). These data demonstrate that the inclusion of space charge may be of great importance in theoretical determinations of efficiency. In Figure (4-2) a plot of efficiency versus current density shows that in the absence of overpotentials we expect behavior similar to that of an ideally Operating fuel cell (see Eq. (4-16)). The ratio. (R). of the voltage loss obtained from Eq. (8-18) to that obtained from integration of Eq. (4-30) obeys the linear relationship R = e + fp. (4-60) The coefficients e and f for Eq. (4-60) are displayed in Table (4-1). Figure (4-3) shows the relationship between R and the current. The increase in R with both p and the current is not surprising since the higher the current or charge_ density the greater the departure from linearity. where higher order solutions are expected to make a greater relative contribution. 56 TABLE 4-1. Least Squares Values of the Empirical Constants Given in the Equations (OV - OIV) = a + bp, n = c + do, and R = e + £0 for Various Values of the Current. Current 2 2 2 3 3 Densiti' a bXIO CX10 dX10 e><10 leO (A-Cm' ) (volts) (volts) 0.00 -0.135 -9.58 0.00 0.000 5.251 10.91 0.05 -0.167 -9.58 3.74 -0.400 5.982 12.40 0.10 -0.201 -9.54 7.19 -0.800 7.043 13.77 0.15 -O.241 -9.58 10.26 -1.169 8.44 16.04 0.20 -0.283 -9.52 13.02 -1.584 10.45 18.31 0.25 -0.330 -9.52 15.30 -1.980 13.01 21.63 0.30 -0.383 -9.52 17.02 -2.388 16.94 25.33 57 Figure (4-2). Dependence of efficiency on current density. (a) Ohm's Law behavior. (b) behavior predicted by Eq. (4-57). 0......0 Amuqv seamen ANIEUIO. 76 Amnmv magmas ANIEoIS 3550 Ecctao end 86 2.0 8d _ _ . _ _ 000 1r .#.on [I II 8.2 I. .II cor: _ _ . _ _ 8 om Kwapwa 77 Either the cathode is cooler than the anode (AT < 0. Q. > 0). or the cathode is warmer than the anode (AT ) 0. 0’ < 0). In either event. because of its nonisothermal properties. the electrolyte prefers to be in the environment of the cathode. Therefore. to produce a current as great as the equivalent isothermal cell requires less ohmic loss across the electrolyte and results in an enhanced efficiency. Exactly the apposite is true for the case of Q‘AT ) 0. Now. regardless of which direction the temperature gradient lies. the nonisothermal properties of the electrolyte dictate that it prefers the environment of the anode. Therefore. to produce the same current as the equivalent isothermal cell requires greater ohmic loss across the electrolyte and results in reduced efficiency. For selected values of Q.AT the efficiency has been calculated using Eq. (5-48). and the results fit to the empirical equation I: - a + le‘AT/T]. (5-53) Although heat of transport data are not available for this electrolyte. the values of the empirical constants a and b displayed in Table (5-1) indicate that the greater the heat of transport of the electrolyte. the greater the potential for enhancing the efficiency. and the more profitable it becomes to exploit the nonisothermal properties of the electrolyte. 78 TABLE 5-1. Least Squares Values of the Empirical Constants Given in n = a + b(Q*AT/T) for Various Values of the Current. .. .. -..M ~——.—-—- ... _ .._.- ——-.—.—.---——-——— - _. ~— _ . ___—_.——————.. I (A—cm'z) __£L__ b (kcal-1 mol) 0.05 3.66 -5.00 0.10 6.98 -10.00 0.15 9.91 1-14.84 0.20 12.42 -20.16 0.25 14.40 -25.00 0.30 15.79 -30.00 CHAPTER 6 NUMERICAL APPROACH A. Introduction The transport equations and boundary conditions described in Chapters 2 and 3 form a set of equations which are. in general. not solvable in closed analytical form. In an effort to verify the accuracy of the pertubation solutions obtained in Chapters 4 and 5 it is reasonable to resort to numerical techniques. Moreover numerical techniques can provide extremely accurate simulation results. The numerical technique of choice is the finite difference method. The details and theory behind the finite difference approach have been, available for some time (Smith [1978] Rosenberg [1969]). The essence of the finite difference approach is the replacement of the continuous space time domain with a discrete space time grid such as that depicted in Figure (641). In doing this the starting set of partial differential equations Eqs. (2-34) to (2-36). is converted to a set of algebraic. or finite difference. equations. In general. the grid spacing is nonuniform both in space and time. The nonuniformity in time allows one to expand the time steps as one approaches the steady state. where dependent variables are essentially constant. and to compress the time spacing at short times. when dependent variables are changing very rapidly. Nonuniformity in the space grid exists for a similar reason. 79 80 Figure (6-1). Nonuniform space and time grid. 81 Ax._, AX O O O O O O O O O O O O O I O 0 O O O O O I O O 0 O O 0 0 e e e 0 e e O O O O O H 1 3+: x _. Figure (6-1) 82 It is expected that the dependent variables will exhibit their greatest spatial dependence in the vicinity of the boundaries. where charged surfaces exist. Therefore it is desirable to increase the density of space points in these regions at the expense of the bulk of the electrolyte. where dependent variables are expected to exhibit very little spatial dependence. B. Finite Difference Formulation Leckey [1981] was the first to transform the correct set of isothermal transport equations to finite difference form. The following discussion derives from that formulation. with appropriate modification for nonisothermal non uni-univalent systems. The dependent variables. denoted by U (U - S. A. E) are specified in space and time with the indices 1 and n respectively. Thus. Ui.n represents the value of U at space column i and time row n with 1 .<. 1 S a. (6-1) H |A 3 IA 0. (6‘2) where R is the number of space points per time row and Q is the number of time rows. There are therefore for J dependent variables. J-R unknowns and JoR equations at each time row. The method for converting Eqs. (2-34) to (2-36) to their algebraic 83 analogs consists of expanding both Ui+l/2.n+1 and Ui-1/2.n+1 in Taylor's series about ui.n+1' The weighted sum of the two series represents the second derivative approximation. and the weighted difference represents the first derivative approximation. with - 2 2 2 (an/3”1,n+1 ' 2 [ U1+1/2.n.+1(“1-1’ + "1,n+1[‘A‘1-1’ ‘ (Ali) 1 “ Ui-l/Z.n+1(A31)2 ]/[(Axi)(Axi_1)(Axi + A‘i-1’l' (6—3) 2 2 (a "/33 )i.n+l ' 4 I “1+1/2.n+1(“1-1’ + Ui-1/2.n+1“‘1’ - Ui'n+1[Ax1_1 — A11] ]/[(Axi_1)(Axi)(Axi_1 + Ax1)]. (a—4) where Ari - x1+1 - xi. and terms of order greater than (Ax)2 are neglected. Similarly. the time derivatives are approximated by expanding Ui.n in a Taylor's series about Ui.n+l and truncating beyond the linear terms. This results in little error if an appropriate time expansion scale is used. The resultant time derivative approximation is where At - tn+l - tn‘ Eqs. (2-34) to (2-36) are converted to finite difference form by repeated application of Eqs. (6-3) to (6-5). The first and last space points require special treatment. Strict 84 application of Eqs. (6-3) and (6-4) to the transport equations requires knowledge of both dependent variables and transport coefficients (ISE. YAE) at the points i - 0. R+1. and transport coefficients (YSS. YAA' YSA' YAS) at grid points i - 1/2. R+1/2. In both cases these grid points lie outside the system boundaries. The treatment of these "image points" requires knowledge of the behavior of the transport coefficients in these regions. These coefficients are in turn dependent upon the concentrations of each species as well as the temperature. Tb determine exactly the nature of the behavior of these coefficients at the boundaries requires that the solution be available. Therefore we are forced to approximate their behavior by noting that for ideal systems the coefficients YSS' YAA’ YSA' and YAS are independent of concentration and can therefore be evaluated at grid points 1 and R instead of at 1/2 and R + 1/2. The electric field at the image points is eliminated through the use of Poisson's equation. The fluxes written in finite difference form are equated to the kinetic boundary flux equations. and both 8 and A at the image points are eliminated. The remaining coefficients YSE and YAE are determined by noting that for an ideal system they are linearly dependent on concentration. In the region of the left boundary we therefore write (OISE/ax)i,l = [ (133) 1,2 - a“) 1,0 ]/2Ax1 s Islsi.2 - Si_o]l2Ax1 + KAlAi-z - Ai-OJIZAxl’ (6-6) where the constants Is and IA are the S and A derivatives of YSE evaluated at the mean system concentrations. From the solutions to si-O 85 and ”i=0 and Eq. (6-6) YSE at the image points is eliminated. The full set of finite difference equations is shown in Appendix E. The finite difference procedure described above has been coded as a Fortran program for the purpose of checking the accuracy of the analytical results obtained in Chapters 4 and 5. A complete listing of the program is included in Appendix F. C. Numerical Results The data available for the doped CeOz electrolyte the steady state voltage loss across the electrolyte has been determined numerically. A comparison of the numerical results versus the analytical results obtained through integration of Eq. (4-30) is shown in Table 6-1. The relative deviation of the two solutions (Adana/A60) increases with the current density. This is not surprising since Ado represents the solution to the linearized problem and as the current density increases the starting equations (Eqs. (2-34) to (2-36)) become more nonlinear. It is for this same reason that the relative contribution of higher order solutions to the potential loss across the electrolyte increase as the current density increases. As we saw in Chapters 3 and 4 the rate constants are dependent upon the behavior of the particular electrode-electrolyte system. In the numerical simulations the bulk potential drop is unaffected by the choice of rate constants. This does 32; mean that the cell performance is independent of the choice of rate constants. since as we saw in 86 TABLE 6-1. Comparison of A¢0 Obtained from Integration of Eq. (4-30) and the Computer Simulation Results for the Doped Solid Electrolyte Ce02-+Ca0.. 1 (A-cm'z) A¢NUM (volts) A¢O (volts) A¢NUM/A¢O 0.00 -0.112 -o.110 1.02 0.05 -0.150 -0.147 1.02 0.10 -0.192 -0.187 1.03 0.15 -O.238 -0.232 1.03 0.20 -0.239 -0.280 1.03 0.25 -O.346 -0.334 1.04 0.30 -O.408 -0.394 1.04 87 Chapter 4 this choice is a reflection of the overpotential loss occurring at the interface. and is the primary cause of power loss. The numerical as well as the analytical results focus only on the bulk behavior which is independent of the properties of the interface. Thus although agreement between the numerical and analytical results based on the transport theory is good they are strictly valid only when no current is flowing. and give increasingly misleading results as the current density increases. The ohmic loss across the electrolyte is defined as the potential drop that would be observed if only pure electric conduction were occurring and the potential drop across the electrolyte were given by name - - I [I/A(S.A)]dx. (15-7) The right hand side of Eq. (6-7) is integrated numerically and a plot of the ohmic potential drop versus current density is shown in Figure (6-2). For a conductor with homogeneous composition the specific conductance of the bulk is defined as (Levine [1978]) 1b = IL/Ad. (6-8) which is got a statement of Ohm/s law but simply the definition of Ab. and applies to all substances. From the plot in Figure (6-2) a value for lb of 0.117 (ohms-cm).1 is extracted by fitting the data to Eq. (6-8). This value compares quite favorably to a value of 0.113 (ohm-cm).1 measured for the same system using alternating current 88 Figure (6-2). Dependence of the ohmic voltage drop on current density. 89 Amlov ouamfib AN! Eol3mco0 “cotmo 0nd 0N0 000 _ _ +- 00.0 -h Im— - ~0— 0nd! 6 1f 8.9. If: 09.0l -r- .4— b — 00.0 (91100) 3110qu 9O techniques (Braunshtein. Tannhauser. and Riess [1981]). The predictions of the transport equations are much more accurate when considering only bulk properties such as lb. An important aspect of the analysis which has been neglected is the effect of overpotential losses on cell performance. Although overpotentials are not experimentally measurable their values can be estimated based upon measurements of the cell impedance 2(0) or admittance 1(0) I 2.1(0) taken over a wide range of frequencies a. Sluyters [1960. 1963. 1964. 1965] has used this technique extensively in his study of aqueous cell polarization phenomena. If one plots the imaginary part of these complex impedances or admittances versus the real part the resulting locus shows distinctive features for certain combinations of circuit components. Thus. these plots are useful for determining the appropriate equivalent circuit for the system under study. Examples of such plots for simple circuits are shown in Figure (6-3). The exact equations for these curves can be derived from a.c. circuit theory (Euler and Dehmelt [1957]). The resistance values in the circuit representations are related to the circular arc intercepts on the real axis. while the capacitance values are related to the frequencies at the peaks of the area. Once the appropriate equivalent circuit for the system has been selected. it is necessary to determine which portions of the equivalent circuit correspond to the bulk electrolyte properties. and which portions of the circuit correspond to the interfacial properties. In their measurements on the daped CeOz platinum paste electrode 91 Figure (6—3). Equivalent circuits and their corresponding admittance plots. The conductance (G) is plotted versus the susceptability (B). 92 AL he 9 in s o R e ft 0 .m 05322: 3 0 no; .................... B .... ( m G o m. R r. + / R C R R / 1 05.3305 3 B 93 electrolyte system. Braunshtein. Tannhauser. and Riess [1981] observed a clearly defined quarter circle. which corresponds to the circuit shown in Figure (6-4) (Bauerle [1969]). The Warburg impedance is the analog for a diffusive process. From the arc intercepts the values of R1 and R2 are obtained. with n, - 3132/(21 + 32) (6-9) no ' 810 (6’10) where R. is the measured resistance for ~19 I and R0 is the measured dynamic resistance for u I 0. The value of R1 thus represents the resistance for current entering the solid electrolyte directly from the gas-electrode-electrolyte interface. The overpotential defined in Eqs. (3-7) and (3-8) is approximated by “L + “a . 131. (6-11) Typical values for R1 and R2 for the Ce02 platinum paste electrode as a function of temperature are shown in Table (6-2). Since R1 >> R2. Eq. (6-9) reduces to Rb . Ru 2 32 (6-12) where Rb is the bulk resistance. It is from these resistance measurements that the literature value of 0.113 (ohm-cm).1 for the specific conductance was extracted. Figure (6-4). 94 Warburg Impedance plot. 95 The "Warburg impedance" is designated by the symbol ~—\A/——- and is equivalent to the infinite R-C line 8 R1 9—; R2 A G /R 45° \/ 1 \\ ,” 1 (1/R1+ 1/R2) Figure (6-4) 96 Inclusion of the overpotentials estimated from Eq. (6-11) results in a far greater loss of power than that predicted by the transport equations alone. A comparison of the power output excluding overpotentials and including them is shown in Table (6-3). and a plot of efficiency versus output current illustrating this same comparison is shown in Figure (6-5). Comparison of these results with those shown in Table (4-2) and Figure (4-2) indicates that neglect of overpotentials is a far more serious error than neglect of space charge in predicting cell performance. 97 Figure (6-5). Dependence of the efficiency on the current density. (+) Neglecting overpotentials. (0) including overpotentials. 98 Aniov seamen ANIEolé 33:00 «cotao 0N0 m. —.0 0 ...0 00.0 00.0 _ 0 _ . _ a m 1+1 ace 0 O O x Li. if. 00.0 x II. 00.0— 00.0, ~_ 001 x K311810043 99 TABLE 6-2- Comparison of the Values of R1 and R2 Measured at High Temperatures for the Doped Solid Electrolyte Ce02-+Ca0. T (°C) R1 (ohms) R2 (ohms) 700 35.71 4.88 750 25.00 3.00 800 18.57 3.27 850 13.21 3.17 900 10.00 4.66 TABLE 6-3. Comparison of Power Densities With and Without Overpotentials. ------ Power (W-cm-2)><102------ _2 Without With I (A-cm ) Overpotential Overpotential 0.00 0.00 0.00 0.05 4.49 3.26 0.10 8.63 3.71 0.15 12.35 1.29 0.20 15.63 -___ 0.30 20.44 .___ CHAPTER 7 THE TIME DEPENDENT PROBLEM A. Introduction Throughout this work we have assumed that the transient period of Operation of the fuel cell (the period from closing the circuit until the steady state is achieved) is short compared to the total operating period of the cell. Tb test the validity of this assumption it is necessary to examine the time dependent behavior of the system within the framework of nonequilibrium thermodynamics. For systems with reversible electrodes the transient period of operation is determined by the thme constant for the establishment of a steady state for mass transport across the bulk electrolyte. The bulk phase is essentially electrically neutral except within a few Debye lengths of the boundaries where kinetic effects dominate. Thus. for estimating the magnitude of the mass transport time constant. we treat the electrolyte as one mutually diffusing species under the influence of temperature and concentration fields. The rate of attainment of the steady state depends upon the magnitudes of the mutual diffusion coefficient. and the thermal conductivity of the solution. We suppress the dependence of this rate upon the thermal conductivity by noting that the thermal relaxation time is generally much shorter than the diffusive relaxation time. Thus we define time zero to be the point at which the 100 101 steady state thermal gradient has established itself but before appreciable chemical diffusion occurs. The transient period of operation of the cell thus becomes a function of the mutual diffusion coefficient. For the doped solid state conductor Ce02 + Ca0. the transport of material occurs via a vacancy mechanism operating in the anion sublattice (Schmalzreid [1974]). We therefore choose as a reference species the Ca+2 dopant which is bound to the lattice positions of the essentially rigid cation framework. which results in a negligible reference velocity. The more nearly perfect the cation sublattice. the closer the reference velocity approaches zero. since there are fewer vacant sites to which the reference species can move. We choose the concentration as the composition variable because it allows us to formulate the time dependent problem directly in terms of the fundamental quantities o and D. and results in an explicit dependence of the composition variable upon these quantities. This in turn leads to estimates of the length of the transient period of cell operation. B. Formulation of the problem Historically the Soret problem has been in a certain sense a zeroth order problem in that electroneutrality has always been assumed. As we saw previously. this assumption can lead to significant errors in the calculation of output voltages when the electrolyte approaches a charged surface such as an electrode. 102 Implicit in the electroneutrality assumption is the requirement that the velocities of the two ions are the same. If the positive ion moved faster than the negative ion. it would leave behind pockets of negative charge while creating pockets of positive charge. and the electroneutrality [assumption would be invalidated. By Eq. (2-1) the faradaio current (IF) becomes. under the electroneutrality assumption. 1F - F E (ziciHvi-vo) .. 12 2 (cizian-vo) = o (7-1) where v m I v1 (i I 1.2) is the velocity of the mutually diffusing electrolyte. The total current defined by Eq. (2-29) then becomes I I s(3E/0t) I e(dE/dt). (7-2) which says that in the steady state the total current must vanish. and that the magnitude and duration of the displacement current is governed by the current produced by the cell as a function of time. Dunlop and Gosting [1959] measured the thermOpower of the thermocell I Pt I Ag I AgN03(aq) | Ag [Pt I. (7-3) and found that after imposition of a temperature gradient across a cell initially uniform in temperature several minutes elapsed before the displacement current vanished. 103 The initial value of the thermopower is taken to be the value of the thermopower measured at the time when the steady temperature field is established. but before appreciable chemical diffusion has occurred. The basis for assuming that the steady temperature field is established before noticeable- diffusion occurs is that the rate of thermal conduction is ordinarily much greater than that of diffusion. The steady state temperature profile is linear since no Joule heating is present. From Eq. (5-49) the entrOpy production for a nonisothermal binary electrolyte system is where the component fluxes are subject to the electroneutrality condition. Ii ' v11”. 1 3 1,2, (7—5) where III is the flux of the mutually diffusing electrolyte. In addition. the driving forces Xi are related through the definition of the chemical potential of a dissociated electrolyte. and X. I -au./ax is the driving force for the mutually diffusing electrolyte. From Eqs. (7-4) to (7-6). 104 e - J.;_ + quq. (7-1) and for the Hittorf diffusion fluxes. -Jn I l-(6un/0x) + l.q(aln T/ax). (7-8) 'Jq I lq.(aun/Ox) + qu(aln T/ax). (7-9) where 1g. - lnq. (7—10) The Onsager coefficients are related to the Soret coefficient. the thermal conductivity. and the mutual diffusion coefficient. the three experimentally accessible quantities required to characterize this type of system. The Hittorf diffusion flux defined by deGroot and coworkers [1953]. is. "I. - (n/covo) [ (z/L) - <1/2) 2 2 + (4/n ) E [cos[(2n+1)uz/L)]/(2n+1) ]exp(hnt) ] (7-34) which corresponds to the solution to the diffusion equation given by Carslaw and Jaeger [1959]. As t 9 I in Eq. (7-54) the steady state solution C 80(x.°) I -(J‘/S-Do)x (7-55) is approached (x I s - (L/2)). From Eq. (4-30) we have. with p = 0. So(x) - I a,oaAo - “103.0 ]x/BA0. (7-56) 113 The solutions given by Eqs. (7-55) and (7-56) are equivalent. This is easily demonstrated by substituting Eq. (4—27) into the right hand side of Eq. (7-56). giving So(x) I (dSo/dx)(xlSn). (7-57) By setting the two solutions given by Eqs. (7-55) and (7-56) equal to one another we regenerate with the help of Eq. (7-57) the diffusion equation —J8 I Do(dSo/dx) (7-58) for an isothermal system. The characteristic time (0) governs the rate at which the steady istate is approached. and is defined by Tyrell [1961] as e = 1.2/(.1200) (7—59) where L is the length of the system and D0 is the zeroth order approximation to the diffusion coefficient. The terms in the sum of Eq. (7-53) converge rapidly and those for which n > 1 can be neglected when tie is sufficiently large. Indeed deGroot [1947] assumed this was always so and neglected higher order terms at all times. If we retain only the first term. Eq. (7-53) simplifies to So(z.t) I 80(z.I) - exp(az+bt)(o11161/0)exp(-t/0). (7-60) 114 Weinberger [1965] defines an upper bound for the error made in approximating an infinite Fourier sum with a partial sum Sm(z). where m is the number of terms included in the sum. For a generalized Fourier series of the form f(x) I a0/2 + E [ ancos(nz) + busin(nz) ] (7-61) this error (lf(z) - Sn(z)|) is given by 1/2 In.) - S.(z)| =[ (1/«1I [f'(z)]2dz - § 1121.1} + 1.1121 ] [ («z/s) — 2 (1111’) II”. (7-62) Application of the above to Eq. (7-60) reveals an upper bound for error of 16 per cent. For purposes of estimating the time required for the system to reach the steady state Eq. (7-60) is sufficient since at t > 50 the error involved in using Eq. (7-60) is less than 1 percent. Eq. (7-60) is. therefore. sufficient for estimating the length of transient behavior in cell operation. It also reveals the effect of the Soret coefficient upon the length of this transient period. From Eq. (7-30) we see that regardless of the algebraic sign of the Soret coefficient or the direction of the temperature gradient. b is always positive in Eq. (7-60). thus the time required to reach the steady state decreases as the magnitude of the nonisothermal effects increases. This decreased warm up time could be economically beneficial provided the nonisothermal effects are large enough. CHAPTER 8 FUTURE WORK To correctly model the time dependent behavior of electrolyte solutions. and gain information concerning single ion properties. it is necessary to know the dependence of the Onsager coefficients upon individual ionic concentrations. This requires making measurements in the first 10"8 s of a liquid junction potential experiment. (Leckey and Horne [1981]). Until such data become available it is difficult to make any statement concerning the validity of any given model in this time domain. Measurements on electrolytes have been carried out on essentially neutral solutions. (Miller [1966]) and necessarily yield information only on the dependence of the Onsager coefficients on the total electrolyte concentration. and forces us to treat the electrolyte as if it were a singly diffusing neutral species. This poses no problem until one reaches a charged interface. where it is necessary to postulate some sort of single ion dependence for the Onsager coefficients. Introduction of a temperature gradient has been suggested as a means of increasing the time available for observing single ion behavior. However this introduces additional parameters (0.1. and KO), raising to six the number of independent observations required to characterize a binary system in this time domain. An additional 115 116 complication is that the temperature dependence of these parameters must also be considered. The introduction of the temperature. and the energy balance equation brings to four the number of dependent variables and linearly independent equations respectively. that are required to characterize this system. The numerical treatment of this problem will require the redemensioning of program "IONFLO" and modeling of the dependence of the above parameters on the temperature and the single ion concentrations. Further experimental data are required in order to more fully develOp this treatment. Alternatively one can treat the electrolyte as one mutually diffusing species. reducing the number of equations and dependent variables to two. This treatment decouples the electric field from the dependent variables. (see Eq. 7-2). thus yielding information only on the temperature and concentration distribution of the system. "IONFLO" can be revised to simulate such a system. and the data necessary for carrying out the simulation. (Miller [1966]. Longsworth [1957]. Snowden and Turner [1960]). are available for liquid electrolytes. Similar data for solid electrolytes are virtually nonexistant. Experiments to determine the temperature and concentration dependence of the diffusion coefficient. Soret coefficient. and thermal conductivity of solids are needed. The electrodelelectrolyte/electrode system is not. in general. symmetrical. The solutions obtained in Chapters 4 and 5 can be further generalized by imposing nonsymmetrical boundary conditions of the type encountered in most real systems. This may result in predictions of 117 significant contributions to the voltage loss from Joule Heating. Another area deserving of more theoretical and experimental attention is the effect of hydrostatic pressure upon transport phenomena in solids. One might expect large effects due to the stress gradients which may be present in many practical situations. APPENDICES and The Y-. q 11 APPENDIX A TRANSPORT COEFFICIENTS in Equations (2-34) to (2-36) are (1/2)[ (011 + D22) ' [(21/22)D12 + (22/21)021] ]L + + + (1/2)[ (011 022) [(21/22)012 (.2/21)0211 ]. (1’2)[ (”11 ”22) ' I‘ll/12’912 ‘ (12/21’9211 ]- ... (1/2)[ (011 022) [(21/22)D12 - (.2121)0211 ][ (F2/3)[ 23122 - z§111 ]. (F2/e)[ 2&111 + 22122112 + 2%122 ]' (1/22122)[ (-Q;/T)(21111 - 22112) + (03/1)(22122 - .1111) ]. <1/22122)[ (QI/T)(z1111 + .2112) + (cg/r)<.1112 + .2122) I. (1-1) 118 119 with K I Nun/Den. (A72) 9 where Num 3 -YSS(zllll + 12112) - YAS(11111 - 12112). (A-3) Den I v1F[ YSS(11111 + 2z122112 + 12122) - YAs(z%122 - 2&111) ]. (Ar4) For electrically neutral regions in an electrolyte. A = 0. and D I [YSSYAE - YSEYASIIYAE' (AIS) with -oS = 05/[r(au/aS)l = [quYAE - YSEYAq1/(DYAE) (Ar6) where D and c are the diffusion coefficient and Soret coefficient of the electrolyte respectively. The a's. B's. and 7's in Equations (4-7) and (4-8) are as ’ [YSAJX ’ YAAJSJIIYSSYAA ‘ YASYSAl' O «A = [YASJX - YsstlllYSSYAA - YASYSA]. 120 as = [ [YSEXAA - YSAYAEJIIYSSYAA - YASYSA] 1(3/22122F)a BA [ [YAEYSS _ YSEYASJIIYSSYAA ’ YASYSA] ](8/22122F)p 75 ‘ [ISAYAq - YSqYAAIIIYSSYAA - YASYSAI' ' ”ASYSq ".YAqusl/ [YssYAA ' YASYSAJ . 0‘ p. I (A-7) APPENDIX B FIRST ORDER SOLUTION FOR THE STEADY STATE ISOTHERMAL PROBLEM Substitution of the zeroth order results into Eq. (4-33) results in the ordinary differential equation (dzfilldxz) - (A3)El = cl: + C2sinh(on) + C3sinh(2on). (3-1) where the constants C1. C2. and C3 are known from the zeroth order solution. In addition to the particular solution which satisfies the right hand side of Eq. (B-l) there is the complementary solution obtained by solving the homogeneous problem. The complementary solution introduces two new constants whose values are fixed by the boundary conditions. From Eq. (4-29). A1 I 0 at the walls. thereby fixing the values of the constants through Poisson's equation with dElldx I O. x = t L/2AA. (B-2) The general first order solution for the electric field is thus given by E1 = T1 + T2 + T3 + T4 (B-3) where T1 = [BAsw/(AOBAOHI [siuh(on)/cosh(Ao§)] - A0: I. 121 122 :3 [aPBA3/(4BAosinh(Ao§))][ (A0.)2sinh(on) - (on)cosh(on) ~ [(Ao§)2 + (Ao§)tanh(Ao§) - 1 ]sinh(A0x) I. r3 - [pE°/(2BAosinh(A0§))][BAA + (psopAs/0A0)1[ (on)cosh(on) - [1 + (AO§)tanh(Ao§)]sinh(on) I. T4 = I(9)2Ao/(6(BA0)281nh2(Ao§))1[BAA + (BAsBso/BA0)1 -[ sinh(2A0x) - 2[cosh(2Ao§)/cosh(Ao§)lsinh(on) I. (3-4) with ° ' [asofiao ’ “aofisol/Bao’ (3'5) E0 - '“Ao/BAO: (B—6) and the Lij' (i.j I S.A) are as defined in Chapter 4. The charge density A1 is obtained from differentiation of Eq. (B-4) in accordance with Poisson's equation. To obtain the first order solution for the total concentration 8. the first order solution for E given by Eqs. (B-3) and (B—4) is substituted into Eq. (4-34). producing (dSIIdx) I Clx + Czsinh(on) + C3sinh(2on) + C4x2sinh(on) 123 + Csxcosh(A0x). (B-7) where the constants C1 to C5 are known from zeroth order results and E1. and are not necessarily the same as the constants given in Eq. (B-l). Integration of Eq. (B-7) produces the first order solution for 8 along with an integration constant whose value is fixed by Eq. (4-28). which states that the total mass of the system is constant. The complete first order solution is Sl=n+T2+T3+T4+T5+T6. (B-8) where T1 3 K1[32 - :2/3le T2 I Kzlxsinh(A0x) - (cosh(on)/Ao) - (cosh(A0§)/Ao) + (2sinh(A0§)/(A3§))l. r3 = x3[ [(on2) - (AO§2)]cosh(on) - [Ssinh(Ao§)/(A3§)] - 3xsinh(A0x) + (Scosh(A0§)/Ao) + (sinh(A0§)tanh(A0§)/Ao) - {tanh(A0§)cosh(on) I. 124 T4 I K4 [cosh(2A0x) - (sinh(2Ao§)/(2Ao§))]. rs = x5 [cosh(on) - (sinh(Ao§)/(A0§))]. rs = x6 I [1 - (Bso/BAO)I[(cosh(on)/Ao) - (s1nh1[0AA + (pAspso/BA011. x2 = [33320.2/(4agpgoco.h(aoc)1. x3 = -[pAss°a/(A30Ao)l[ (E°/2)[Bss — (BAsBsolflAo)1 - (BAso/BA0)I. K. = [BAsa(E°)2/(A30Ao)l[ 0.. + ipss - (BAsBso/PAo)ll2/61 + [BAsBso/BAOJ[tanh(Ao§)/(Ao§)]. (3-14) and 32(. = 0) a K1(on)2sinh(on) + K2(on)cosh(on) + K3sinh(on) + [4,3 + x53, (B-15) where. moreover. x1 = [BAsa/ZBA0]2[BSOE9/(Agcosh(AO§))1. K2 2 [pAsaE°/(A30A0cosh(aoc))][ (aBss) + (BSOBAAEOIZ) - [(BAsBso/ZBAO)2(aBSOE°)] I. 127 K. . 1as°0As/(Agco.h(aoc))1[ tapsoaAs/4niolt1 + (Aoc)2 + (Ao§)coth(A0§)] - (EOaso/0A0113AA + (BAsBso/BAo>1 -[1 + (Act/2)coth(A0§)] + [BSA + (BSOBSS/BA0)]E° - [afiss/BAol I + [2taE°9sofias/(A6flaosinh1A05)>1 [ (s°/2)tass - (flASBso/BA0)] - (oBAs/BAO) I. K4 3 (aE°/3)[Bss - (BsoflAs/BAO)][ (E°/2)[Bss _ (BsoBAs/BA0)1 ‘ (“Baa/3A0) 1' K. = -[aE°/A%][Bss — (BAsBso/BAO)][ [(Ao§)2/6][Bss - (BASBSO/BA0)] + [BASBSOIIAOBA0§)lt‘nh(A0§) ] + [BAss°a/A31[ E°(Bso/5Ao)[BAA - (ass - (BAsBso/BA0)) - BSA] + (zaaAsasoiafio) ]. (3-16) The preceding results represent the second order contributions to the electric field and total concentration due to bulk behavior only. Integration of Eq. (B-l3) from x I -§ to +§ produces A‘2(p'0) = (2E°§3/3)(oBAsIBAo)2. (3-17) 128 where terms of order 62 and below have been neglected. It is from Eq. (B-17) that the voltage drops corresponding to the A/SIII I 0 entries in Table (4-1) are calculated. Inclusion of terms linear in p enables us to examine .the effect of non electrically neutral behavior upon the voltage loss across the electrolyte. and results in the expression Adz - A62(pI0) - (op/BAO)(BAs§/BA0)2[ (3o/8) - (E°Bso/3) I, (3.13) where terms of order 52 and below have been neglected. The second order contributions to the voltage drop as a function of charge density in Table (4-1) are computed from Eq. (B-18). APPENDIX C DIMENSIONLESS VARIABLES USED FOR COMPUTER SIMULATION Dimensionless variables are denoted by an overline. The following list contains dimensionless variables used in addition to those introduced in Eq. (4-16). The definition of each variable is given in the comment section of the program ”IONFLO." cpl I 8(3I0.KIO). EPS I spm/so. .t- . tIYnAEmll CON ‘ 80/(22122F)e GDE a CpflYmAE/(22122F), 20 =- -Q(C0N)/[eo(M)sm], 21 = -I(c0N)/[(AA)sngEepml. E3 = -Ae(C0N)/[2(AA)ZSn10-7]. 129 130 213 - Yij(CON)/[(AA)ZGDE]. 1.1 = S.A.q if? - kij(CON)/[(AA)GDE]. i = S.A; j = L.0.A.R. APPENDIX D JOULE HEATING Neglect of the Joule heating term in the energy balance equation (Eq. (5-10)) resulted in a linear temperature gradient. The resultant steady state equations for nonisothermal transport (Eqs. (5-24) and (5-25)) thus assumed a form identical to their isothermal analogs (Eqs. (4-26) and (4-27)). Inclusion of the Joule heating term in the energy balance equation transforms Eq. (5-24) to (don/dx) - (A3)Eo = -(l/EPS)[ u'Ao - (21on/(10Korn)) -[IL/2§]2 I, (0—1) which reduces to Eq. (5-24) when I I 0. The general solution to Eq. (D-l) is £0 a (IL/2:12127on/(xoxormaAo)1 - (aAo/0A01 + Clcosh(A0x) + Czsinh(on). (D-2) where the constants C1 and C2 are determined by the boundary conditions given by Eq. (4-29). Under the conditions of Eq. (4-32). which states that the charge density is equal but apposite in sign at either end of the cell. the constants become 131 132 C. = p/[A0(EPS)sinh(Ao§)]. (0-3) C2 ’- 27A0(IL)2/ I (2§)210K0TnAoflAocOSh(A0§)] e (D-4) With C1 and C2 substituted into Eq. (D-2) the zeroth order expression for E0 _reduces to Eq. (4-30) when the current density vanishes. Integration of Eq. (D-2) across the electrolyte produces the same expression for Ado as we obtain from Eq. (4-30) regardless of the value of the current density. This does not mean that the voltage loss across the electrolyte is independent of the current density: rather. it means that for the purpose of calculating the voltage loss across the electrolyte it is not necessary to include the Joule heating term in the energy balance equation. provided that the charge density distribution is symmetrical. If Eq. (4-32) is not valid (pL # -pR), then the constants C1 and C2 defined by the boundary conditions will change. thus altering the final expression for_Aeo obtained by integration of Eq. (D-2). We would therefore expect Joule heating to be important in the determination of voltage loss for systems employing different kinds of electrodes where nonsymmetrical charge distributions are expected. APPENDIX E TRANSPORT EQUATIONS IN FINITE DIFFERENCE FORM The following definitions are useful in simplifying the finite difference form of the transport equations. A - 1/(Axi + Axi-1). P - 4/(Axi + Axi-1)2. Y0v ' YUV(si.n+1' Ai.n+l)' Yfiv ' YUV(Si+l/2.n+l' Ai+1/2.n+1’» Yfiv ' Yuv(31-1/2,n+1: Ai-l/2.n+l)' YUE ‘ YUE(Si+l.n+l' Ai+l.n+l)' YUE ' YUE(si-l.n+l' Ai-l.n+l)' YUE ' YUE(si'n+1' Ai.n+1)' Rat I Ali/Axi-1. Sub I Ali - Ali-1. Amul I (Axi)(Axi-1). K I e/(ZzlzzF). (E-l) where (U = S.A and V = S.A.q). with Ui.j referring to the value of U at space point i and time row j. Equation (2-34) becomes for 2 S i i er. .. + .. Si,n/At 3 -PYSSSI‘1,n+1 + [P(YSS + Yss) + l/Atlsi'n+1 133 134 _ - - + PY§ssi+1.n+1 ' PYsaAi-1,n+1 + P‘YSA + YSA)Ai.n+1 + + [K(Sub)/Amul]YSEEi.n+1 - pyqui+1'n+1 _ PYSqTi-1 + P(Y§q + qu)Ti-n+1' (E—Z) Equation (2-35) becomes. for 2 S i g 321, Ai,n/At = - PYKssi-1,n+1 + p(y;s + YAS)si.n+1 _ - — + pYXssi+1,n+1 ' PYAAAi-l.n+l + (PYAA + PYAA + 1’At)Ai,n+1 + ‘ + (AK)(R.t)YAEEi+1.n+1 ‘ (AK/R1I)Ei_1.n+1 + [K(Sub)/Amul]YAEEi'n+1 ’ PYXqu-1,n+1 + P(YAq - + YAq’Ti.n+1 ‘ PYAqTi+1,n+1- (E—3) Equation (2-36) becomes. for 2 S i S Rel. (I/e) + Ei.n/At = Ei'n+1[(1/At) + YAE] - (AYAs/K) °[(R‘t)si+1.n+1 ' (31-1'n+1/R8t)1 ‘ (AYAA/K)[(R‘t)Ai+1.n+l - (Ai_1'n+1/Rat)] - AYAq[(Rat)Ti+1,n+1 - (Ti-1,n+1/Rat)l. (E—4) 135 The rate constants enter through the elimination of the values of S and A at the image points i I 0. 3+1. The assumption of a linear temperature gradient being established at time zero. with T1 a T0 and Tk+1 I Th. results in the cancellation of all temperature terms from the transport equations at the boundaries. Thus for i I 1. Eq. (2-34) becomes si'n/At + 4A[ (KSL)SL + (KAL)AL] = -2PYSSSZ.n+1 _. 2PYSAA2,n+1 + [(1/31) + 2PYSS + Anursonsl,n+1 + [szSA + 4A(aso) + + YSEIA1.n+1 + 2A“31331.1“: + 2A1‘YSIE.1‘31.11+1- (5'5) For i I 1. Eq. (2-35) becomes Ai'n/At + 4A[(KSL)AL + (KAL)SL] = -2PYASSZ'n+1 - 2PYAAA2’n+1 + [ZPYAS + 4A(KDO)]sl,n+1 + [(llAt) + ZPYAA + 4A(KSO) + YAE] + (A1.n+1) + ZAKYAEEI,n+1 + ZAKYAEE1,n+1° (8’6) For i I 1. Eq. (2-36) becomes (I/a) + (31.n/At) - (1/K)[(KAL)SL + (KSL)AL] = (El'n+1/At) - (1/K)[(KDO)Sl'n+1 + (KSO)A1'n+1]o (E'?) 136 For i = R. Eq. (2-34) becomes 81.n/At + 4A[(KSR)SR + (KAR)AR] = ’ 2PYSSSRPI,n+1 - ZPYSAAR?1,n+1 + [zPYSA + “(mm + YSEMRm-t-l + [(l/At) + 2PYSS + 4A(KSA)] 'SR.n+1 ' ZAKYszEk.n+1 ' ZAKYEEER,n+1o (E-B) For i = R, Eq. (2-35) becomes A1.n/At + 4A[(KAR)SR + (KSR)AR] 3 -2PYASSR71,n+1 - ZPYAAAR-1,n+1 + [(l/At) + ZPYAA + 4A(XSA) + YAEJAR,n+1 + [szA8 + 4A(KAA)] + 08R.n+1 - 2AKYAEER,B+1 ‘ ZAKYAEER,n+1° (8’9) For i = R. Eq. (2-36) becomes (Us) + (Emu/At) + (unuxmsn + (xsmm = (Emma/At) + (1/K)I(KAA)SR,H+1 + (KSA)AR.n+1] (E-IO) nnnnnnnnnnnnnnnnnnnnnnnnnnnnnonnnnnnnnnnn APPENDIX F LISTING OF PROGRAM “IONFLO” The following is a listing of the program ”IONFLO”: PROGRAM IONFLO (INPUT,OUTPUT,TAPE 6I'OUTPUT,TAPE 50-INPUT) PROGRAM IONFLO SOLVES THE SYSTEM OF EQUATIONS, GENERALIZED NERNST- PLANCK, POISSON, DISPLACEMENT CURRENT, FOR THE TRANSPORT OF A BINARY ELECTROLYTE IN A FINITE SPACE DOMAIN. DERIVATIVES ARE APPROXIMATED USING FINITE DIFFERENCES. 2 TYPES OF BOUNDARY VALUE PROBLEMS ARE POSSIBLE. KEY-I IS FOR THE DIFFUSION POTENTIAL KEY-2 IS FOR KINETIC ELECTRODE MODELING BOTH EMPLOY A NON-UNIFORM SPACE MESH WHICH DEPENDS ON THE MEAN CONCENTRATION OF THE ELECTROLYTE AND THE INPUT PARAMETER RANGE. Q INPUT ALL INPUT IS READ THROUGH SUBROUTINE AIO. THIS SUBROUTINE SHOULD BE INSPECTED TO DETERMINE PROPER ORDER AND FORMAT OF INPUT PARAMETERS. THE FOLLOWING IS A LIST OF INPUT VARIABLES AND THEIR MEANING. TITLE DESCRIPTION OF SIMULATION KEY I FOR DIFFUSION POTENTIAL SIMULATION 2 FOR KINETIC ELECTRODE MODELING R NUMBER or SPACE POINTS UP T0 100 Q NUMBER or TIME ROWS ISET NUMBER OF EQUAL INCREMENT SPACE POINTS AT EACH ELECTRODE OR SIDE OF INTERFACE IJK FRACTION OF TIME ROWS TO BE WRITTEN ISTOP NUMBER OF TIME ROWS TIME STEP SIZES ARE CONSTANT MMM SPECIFIES A TIME ROW FOR WHICH MANY INTERNAL PARAMETERS WILL OUTPUT AND SHOULD ONLY BE USED FOR DEBUGGING. WHEN NOT IN USE SET MMM GREATER THAN Q. NPLT FRACTION OF TIME ROWS THAT ARE WRITTEN THAT ARE TO BE PLOTTED 137 nnnnnnnnnnnnnnnnnonnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn IOFF ISTOP RANGE DT TOL ZP ZM ZQ ALF A GDE SIZE 138 I FOR COMPLETE ITERATION PROCEDURE 0 FOR NO ITERATION O WRITES OUTPUT ONLY FOR LAST ITERATION OF EACH TIME ROW SPECIFIED BY IJK I WRITES OUTPUT FOR EVERY ITERATION OF EACH TIME ROW SPECIFIED BY IJK LAST TIME ROW THAT CURRENT IS TURNED ON FOR DETERMINES RATE AT WHICH GRID SPACING EXPANDS FROM SMALL STEP SIZES TO LARGE ONES. INITIAL TIME STEP SIZE IN SECONDS DETERMINES RATE OF TIME STEP SIZE INCREASE. MAXIMUM ALLOWABLE FRACTIONAL CHANGE FROM ONE ITERATION TO THE NEXT FOR PROGRAM TO PROCEED TO THE NEW TIME ROW CHARGE NUMBER FOR POSITIVE ION CHARGE NUMBER FOR NEGATIVE ION FINAL CHARGE ON ELECTRODE IN C/CM**2 TIME CONSTANT FOR CHARGE INJECTION CHARACTERISTIC LENGTH IN CM MEAN ZDE VALUE IN MOL*S*C/G*CM**3 LENGTH OF SYSTEM IN CM EPM PERMITTIVITY OF A VACUUM IN C*C*S*S/G*CM**3. SM MEAN ELECTROLYTE CONCENTRATION IN MOLES/CM**3 KDO,KSO.KDA,KSA. RATE CONSTANTS FOR TRANSPORT ACROSS SR’SL, C(I) ELECTRODE-ELECTROLYTE INTERFACES. DR.DL. RESERVOIR CONCENTRATIONS. THROUGH C(lO) EXPANSION COEFFICIENTS FOR ONSAGER COEFFICIENTS AND ACTIVITY COEFFICIENT DERIVATIVES. C(I) AND c(2) FOR L++, c(3) AND C(A) FOR L--, C(5) ANO C(6) FOR L+-, C(7) ANO C(8) FOR M, C(9) ANO C(Io) FOR 3 UNITS ARE S*MOL*MOL/(G*CM**3)-FOR ONSAGER COEFFICIENTS. M AND a ARE IN CM**3/MOLE. EXPANSIONS ARE OF THE FORM: L++ - CI*CP + C2*CP**I.5 L-- - C3*CM + CA*CM**I.5 L+— - C5*S**l.5 + c6*s**2 nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnonnnnnnnnnnnnnnn 139 M - C7*S - C8/S**(0.5) B - C9 + ClO*S NAME(I,2,3) TITLES FOR PLOTS JNAME(I,2.3) TITLES FOR PLOTS POINT(I,2,3) SYMBOLS FOR PLOTS OUTPUT ALL OF THE INPUT PARAMETERS. THE OEBYE LENGTH OF THE SYSTEM, AND THE TRANSPORT COEFFICIENTS IN OIMENSIONLESS FORM ARE WRITTEN BEFORE THE FIRST TIME ROW IS EXECUTED s SUM CONCENTRATION. WRITTEN FOR EVERY TIME ROW AS INOICATEO BY IJK AND EVERY SPACE POINT AS INDICATED BY IW. UNITS ARE MOLES/CM**3 . . SM IS SUBTRACTEO FROM S FOR THE OUTOUT D CHARGE DENSITY LABELEO RHO ON OUTPUT. WRITTEN FOR EVERY TIME ROW AS INOICATEO BY IJK AND EVERY SPACE POINT As INOICATEO BY IW. UNITS ARE C/CM**3 E SAME COMMENT As FOR 0. UNITS ARE VOLTS/CM PSI SAME COMMENT AS FOR D. UNITS ARE VOLTS RELATIVE TO THE POTENTIAL AT X-O. CP POSITIVE ION CONCENTRATION CM NEGATIVE ION CONCENTRATION x SPACE POSITION IN CM RELATIVE TO NEAREST WALL. T TIME IN UNITS OF SECONOS EJ CELL POTENTIAL DIFFERENCE IN VOLTS ZI DIFFERENCE IN CURRENT DENSITY BETWEEN ARBITRARY TIME T AND STEADY STATE, AMPS/CM**2. CLE TOTAL CHARGE IN SYSTEM IN COULOMBS RCOND SEE LINPACK DOCUMENTATION IF SUBROUTINE SGBFA IS USED BY IONFLO RCOND IS MEANINGLESS RTEST ARBITRARY CUTOFF VALUE TO TELL COMPUTER THAT MATRIX IS NEARLY SINGULAR AND TO STOP ATTEMPTING A SOLUTION. A VALUE OF IO**(-25) HAS BEEN USED SUCCESSFULLY. EPS RATIO OF PERMITTIVITY OF ELECTROLYTE TO THAT OF A VACUUM. SO INITIAL CONDITION ON S. DO INITIAL CONDITION ON D. nnnn nnnnnnn nn nnnn no 15 13 no 140 ED INITIAL CONDITION ON E. QDELT IS THE QUANTITY Q*(TR - TL)/TM IN CALORIES PER MOLE. VI AND v2 ARE THE STOICHIOMETRIC COEFFICIENTS OF SPECIES I AND 2. DELT IS THE TEMPERATURE DIFFERENCE, TR-TL, IN KELVIN. INTEGER R.RR.Q DOUBLE XER.DEL.SIZE.DEBL DIMENSION B(3OO),IPVT(3OO),BB(3OO),2(3OO) DIMENSION Y(300),DEL(IOO),XER(IOO).YY(300) DIMENSION ABD(I6.300) DIMENSION NPT(IO),JNAME(IO).NAME(IO) COMMON /ELEC/ EPS,ZI COMMON /GREEK/ c(2h) COMMON /PARAM/ KEY,R,Q,DT,F,TOL,RANGE,ISET COMMON/UNITS/SM,DF,GDE,EPM,A,ZP,ZM,CON,TM,DIM,CREF,ZREF COMMON /BOX/ IW,NPLT,IJK COMMON/ICOND/SO,DO,EO COMMON /PLTPTS/POINT(IO) COMMON /CHRG/ ZQ,ALF COMMON/MATRIX/RTEST,BCOND SUBROUTINE AIO READS THE INPUT DATA AND CONVERTS VARIABLES T0 DIMENSIONLESS FORM. MOST OF THE INPUT DATA IS ALSO WRITTEN CALL AIO (IOFF,ISTOP,SIZE,DEBL,NPT,NAME,JNAME.MMM) IF BCOND IS ZERO THEN THE PROGRAM WILL SOLVE ONLY THE STEADY STATE PROBLEM. IF (BCOND .EQ. O.O)GO TO 335 SUBROUTINE GRIDXX, WHERE XX IS PR FOR ONE BLOCKING ELECTRODE, RR FOR TWO BLOCKING ELECTRODES, AND RR FOR LIQUID JUNCTION SIMULATIONS, GENERATES THE NONUNIFORM SPACE MESH ACCORDING TO R. GAUGE, AND KEY. DEL IS THE INCREMENT VECTOR AND XER HOLDS THE SPACE POSITIONS OF THE GRID PTS. CALL GRIDPR (DEL,XER,SIZE.DEBL) SHIFTING ORIGIN FOR KEY-l PROBLEM IF (KEY .EQ. 2) GO TO I3 DO I5 I-I,R XER(I)-XER(I)-(O.5*SIZE) CONTINUE CONTINUE RR-3*R TIME-0.0 RCOND-I.O ZI - 0.0 M - o KOUNT IS THE NUMBER OF ITERATIONS PERFORMED AT EACH TIME ROW. KOUNT-O n n on On nnn nnnnn an 18 nnnn 21 nnnnn nnnnnnn 141 SUBROUTINE IC ASSEMBLES THE INITIAL CONDITIONS FOR 5.0, AND E. CALL IC (KEY,Y,R) OUTPUTTING THE INITIAL CONDITIONS CALL lANDW(Y,TIME,KOUNT,RCOND,XER,NPT,NAME,JNAME,M,KFLAG,SIZE) SAVING F AS FF FOR USE AFTER CHARGING IS COMPLETE. FF - F BEGINNING OF TIME LOOP DO I7 M-I,Q ZI - ALF*ZQ*EXP(-ALF*TIME) IN GENERAL TIME(N.F,DT)'DT*((I+F)**N-I)/F TIME-TIME+DT SUBROUTINE RHSVO COMPUTES THE PORTION OF THE RHS-VECTOR GENERATED BY THE SOEUTIONS FROM THE PREVIOUS TIME ROW. RHSVO IS CALLED ONLY ONCE FOR EACH TIME ROW. THE RHS-VECTOR IS RETURNED AS B. CALL RHSVO (R,OEL,DT,Y,B,M,MMM) THE UNCHANGING PORTION OF B IS SAVED AS BB FOR USE WITH EACH SUBSEQUENT ITERATION. DD 18 I-I,RR BB(I)-B(I) CONTINUE SUBROUTINE COEMAT ASSEMBLES THE COEFFICIENT MATRIX, ABD, IN BAND STORAGE FORM. CALL COEMAT (R,DEL.DT,Y,ABD,M,MMM) IF RCOND LESS THAN RTEST THEN MATRIX WILL BE NEAR SINGULAR. RTEST IS SUPPLIED BY THE USER BY OBTAINING THE CONDITION OF THE MATRIX NEAR THE STEADY STATE. IF (ABS(RCOND) .LE. RTEST) GO TO 333 SUBROUTINE RHSVN COMPUTES THE PORTION OF THE RHS-VECTOR DUE TO THE SOLUTION TO THE PREVIOUS ITERATION. THIS VECTOR IS RETURNED AS B. THE VECTOR B+BB IS THE COMPLETE RHS-VECTOR FOR ANY GIVEN ITERATION. CALL RHSVN (R,DEL,DT,Y,B) DO 19 J-I,RR B(J)-B(J) + BB(J) nnnnnn nnnn nnnnnn 20 nnnnnnn anon noon 27 nnnn 142 YYIJ)-Y(J) CONTINUE LDA IS THE LEADING DIMENSION OF ABD. ML IS THE NUMBER OF DIAGONALS ABOVE THE MAIN, MU IS THE NUMBER OF DIAGONALS BELOW THE MAIN. LDA-l6 ML-S Mu-S SGBCO FACTORIZES THE MATRIX ABD. RCOND IS RETURNED AS I/COND. WHERE COND IS THE CONDITION OF THE MATRIX. IF RCOND IS NOT NEEDED. SUBROUTINE SGBFA MAY BE USED TO DECREASE EXECUTION TIME JOB-O CALL SGBCO (ABD,LDA,RR,ML,MU,IPVT,RCOND,Z) SGBSL SOLVES THE MATRIX EQUATION ABD*Y-B. WHERE ABD IS THE BAND MATRIX FACTORIZED BY SGBCO, B IS FORMED FROM RHSVO AND RHSVN. AND Y IS THE SOLUTION VECTOR RETURNED AS B. THE ENTERING RHS-VECTOR. B. IS DESTROYED. CALL SGBSL (ABD,LDA,RR.ML,MU,IPVT,B,JOB) DO 20 K-I,RR Y(K)-B(K) CONTINUE CMPR COMPARES THE CURRENT ITERATION SOLUTIONS WITH THE PREVIOUS ONES. IF THE RELATIVE CHANGE IS LESS THAN TOL FOR EACH VARIABLE. KFLAG IS RETURNED AS ZERO AND THE MAIN PROGRAM CONTINUES TO THE NEXT TIME ROW. OTHERWISE KFLAG-I AND ANOTHER ITERATION IS PERFORMED. CALL CMPR (Y,YY,TOL,R,KFLAG) THE NEXT CARD REDUCES EXECUTION TIME FOR SMALL PERTURBATIONS BY ELIMINATING THE ITERATION PROCEDURE KFLAG-IOFF*KFLAG KOUNT-KOUNT+I IF KOUNT REACHES II ITERATION IS PROBABLY DIVERGING SO PROGRAM IS TERMINATED. IF (KOUNT .EQ. II) GO TO I6 DO 27 J-I,RR Y(J)-YY(J) CONTINUE IF M IS NOT EQUAL TO I OR A MULTIPLE OF IJK. THE OUTPUT AT THAT TIME ROW WILL NOT BE WRITTEN nnnnn nnnn nnnnnnn nnnnn II IA I2 I7 16 I0 333 33“ 335 143 L-MOD (M.IJK) IF (L .NE. D .AND. M .NE. I) GO TO IA IF (IN .EQ. 0) GO TO II IANDW CALCULATES THE CHARGE AND MASS IN THE SYSTEM, THE POTENTIAL AS A FUNCTION OF X, WRITES EVERY IJK TH TIME ROW AND PLOTS S.D.AND E AS A FUNCTION OF X. CALL IANDW(Y,TIME,KOUNT,RCOND,XER,NPT,NAME,JNAME,M,KFLAG,SIZE) CONTINUE CONTINUE IF (KFLAG .EQ. I) GO TO 2I IF (L .NE. 0 .AND. M .NE. I) GO TO I2 IF (IN .NE. 0 ) GO TO I2 CALL IANDW(Y.TIME.KOUNT,RCOND,XER,NPT,NAME,JNAME.M,KFLAG,SIZE) CONTINUE DT IS INCREMENTED AT EACH TIME ROW ACCORDING TO THE VALUE OF F THAT IS SPECIFIED IN THE INPUT DATA F KEEPS TIME STEP SIZES CONSTANT WHILE ELECTRODE IS BEING CHARGED IF (M .LE. ISTOP) F's 0.0 AFTER ELECTRODE CHARGE IS COMPLETED TIME STEP SIZES ARE GRADUALLY INCREASED IF F IS SET POSITIVE IF (M .GT. ISTOP) F - ABS(FF) IF F IS SET BETWEEN 0 AND -I TIME STEP SIZES WILL BE DECREASED BY ABSIFF) FOR TIME ROW ISTOP+I AND THEN GRADUALLY INCREASED THEREAFTER AS USUAL IF (M .EQ. ISTOP .AND. FF .LT. O.O) DT - ABS(FF)*DT DT-(F+l.0)*DT KOUNT-0 CONTINUE STOP CONTINUE WRITE (6I,IO) FORMAT (Ix.28HSOLUTION VECTOR Is DIVERGING) GO To 335 NRITE (61.33h) FORMAT (Ix,23HMATRIx IS NEAR SINGULAR) CONTINUE END SUBROUTINE AIO (IOFF,ISTOP,SIZE,DEBL,NPT,NAME,JNAME,MMM) nnnnnn 144 ’SUBROUTINE AIO READS IN ALL THE DATA NEEDED BY IONFLO. ALL INPUT DATA IS ALSO WRITTEN BY AIO. ALL INPUT FORMAT IS EITHER I5 OR EIO.3 EXCEPT NAME AND JNAME WHICH ARE A5 AND TITLEI AND TITLEZ WHICH ARE AIO. DIMENSION NPT(IO),JNAME(IO),NAME(IO) REAL II,I2,ITOT DOUBLE SIZE.DEBL REAL KSL,KDL,KSO,KDO.KSA.KDA.KSR.KDR INTEGER R,Q COMMON/CURRENT/II,I2 COMMON /ELEC/ EPS,2I COMMON /GREEK/ C(2A) COMMON/UNlTS/SM,DF,GDE,EPM,A,ZP,ZM,CON,TM,DIM,CREF,ZREF COMMON /PARAM/ KEY,R,Q,DT,F,TOL.RANGE,ISET COMMON /PLTPTS/ POINT(IO) COMMON /BOX/ IW,NPLT,|JK COMMON/ICOND/SD,DO,EO COMMON /CHRG/ zq.ALF COMMON/RATES/KSL,KDL.KSO,KDO.KSA,KDA,KSR,KDR COMMON/BNDRY/SL.SR,DL,DR COMMON/MATRIX/RTEST.BCOND COMMON/DIFF/DI.Dz,OELTA COMMON/THERM/ODELT,VI,V2 DOO,AND DA ARE THE DIFFERENCE CONCENTRATIONS AT x--L/2 AND x-+L/2 RESPECTIVELY READ (50,37) TITLEI,TITLE2 READ (50,I0) KEY,R,Q,ISET,IJK,NPLT,IOFF,IW,MMM,ISTOP,RANGE,DT,F READ (50.13) ZP.ZM.ZQ,ALF,EPS,TOL,TM,RTEST READ (50,13) EPM,zREF,SIzE,CREF,SO,DO,EO,BCOND READ (50,I3) CRESL2,CRESR2.DL.DR READ (50,l3) II,ITDT,DI.Dz,DELTA.DOO.DA READ (50,13) (C(I),I-I,IO) READ (50,13) QDELT,VI.V2 READ (50,II) (POINT(J),J-I,3) READ (50,12) (NAME(J),J-I,3) READ (50.12) (JNAME(J).J-I.8) CRESLI - (2*ZP*ZM*DL - CRESL2*ZM - CREF*ZREF)/ZP CRESRI - (2*ZP*ZM*DR - CRESR2*ZM - CREF*ZREF)/ZP SL - O.5*((CRESL2/ZP) - (CRESLI/2M)) SR - O.5*((CRESR2/ZP) - (CRESRI/2M)) SM - (SL+SR)/2.0 SL - (SL-SM)/SM SR - (SR-SMI/SM DR - DR/SM DL - DL/SM DOO - DOD/SM DA - DA/SM DF--ZREF*CREF/(ZP*ZM*2.0) FIII-C(I)*(ZM*(DF-SM)) FIIZ'C(2)*(((ZM*(DF-SM)))**I.5) FIZI-C(5)*((((ZP+ZM)/2.0)*DF + ((ZP-ZM)/2.0)*SM)**I.5) \ 145 Flzz-C(6)*((((ZP+ZM)/2.0)*DF + ((ZP-ZMI/2.0)*SM)**2.0) F22]-C(3)*(ZP*(DF+SM)) F222-C(h)*(((ZP*(DF+SM)))**l.5) GII-(FIII + Fll2)*ZP*ZP GI2-2.0*ZP*ZM*(FIZI + FI22) GZZ-ZM*ZM*(F221 + F222) GDE-(96h85.0/(2.0*ZP*ZM))*(GII + GI2 + G22) SC - (ZP*ZM*SM/2.0)*(((ZP+ZM)*(DF/SM)) + (ZM-ZP) + + ((ZREF*ZREF)*CREF/(ZP*ZM*SM))) DEBL-((8.3)hE7*EPS*EPM*TM)/(2.0*(96h85.0**2.0)*SC))**0.5 IO II 12 13 In 15 I6 22 25 26 27 28 A-DEBL /5.0 KSL- (DZ+DI)/(2.0*DELTA) KDL - (DZ-DI)/(2.0*DELTA) KSA I KSO KSR KDA KDO KDR I WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE FORMAT FORMAT FORMAT FORMAT FORMAT FORMAT FORMAT FORMAT FORMAT FORMAT FORMAT FORMAT KSL KSL KSL KDL KDL KDL (61.33) (6I.28) (61.35) (6I,32) (61.31) (6I.ho) (6I,AI) TITLEI,TITLE2 R,Q,KEY DT,F,TOL,RANGE,RTEST A.ALF IJK,ISET,IW,IOFF,ISTOP,NPLT GDE,EPM,SM,DF.ZQ,SIZE CREF.EPS (6I,38)KSO,KDO,KSA,KDA (6I,83)KSL,KDL,KSR,KDR (6I.53)DELTA (61.27) (61.26) (61.25) (6I.IA) (61.15) (6I,I6) (61.29) C(I).C(2) C(3).C(II) C(5).C(6) C(7).c(8) C(9).C(IO) DEBL.TM (6I.68)SL.SR.DL.DR (IOI5,3EIO.3) (3A1) (8AIo) (8EIO.3) (///.2x.I9H M AND B EXPANSIONS./) (SXQBHM -sE]o°395H/S T 9E100397H*S**]°59/) (3X,3HB -,El0.3,3H + ,EIO.3.2H*S,/) (//.2X,31H ONSAGER COEFFICIENT EXPANSIONS,/) (3X,5HLPM -,EIO.3,IZH*(S)**I.5 + ,El0.3,7H*(S)**2./) (3X,5HLMM -,EIO.3,IIH*(S + D) + ,EIO.3,I3H*(S + D)**l.5,/) (3X,5HLPP -.EIO.3,IIH*(S - D) + ,EIO.3,I3H*(S - D)**I.5./) (/,IX,|A,IOH SPACE PTS,AX,IA,IOH TIME ROWS, I2x.6H KEY -.I2.//) 29 FORMAT (////,2x,I5H DEBYE LENGTH -,EIO.3,7x,5H TM -,EIO.3,//) 3I FORMAT (2x,6H IJK -,I3,8x,7H ISET -,I3,6x,5H IW -,I3,on,7H IOFF . I,I3,on,8H ISTOP -,I3,on,7H NPLT -,I3,//) noon 0 32 33 35 53 68 II9 37 38 83 AD AI 50 146 FORMAT (2x,AH A -,EIO.3,3x,7H ALF =,EIO.3,//) FORMAT (IHI.6OX.2AIO.///) FORMAT (2x,AH DT-,EIO.3,3x,3H F=,EIO.3,3x,5H TOL=,EIO.3,3x,7H RANG IE-,EIO.3,3X,7H RTEST-,EIO.3,//) FORMAT (3X,7HDELTA -,EIO.3,//) FORMAT (//,3X,AHSL -,EIO.3,3X,AHSR =,EIO.3,3X,AHDL =,EIO.3, + 3x,hHDR -.EIO-3.//) FORMAT (IX,6HCOMP -,EIO.3,3X,5HSUM =,EIO.3,//) FORMAT (ZAIO) FORMAT (3X,5HKSO '.EIO.3,2X,SHKDO =,EIO.3,2X,5HKSA =,EIO.3,. + 2X,5HKDA t,EIO.3,/) FORMAT (3X,5HKSL -,EIO.3,2X,5HKDL =,EIO.3,2X,5HKSR =,EIO.3, + 2X,5HKDR I'=.EIO.3,/) FORMAT (2X,5H GDE-,EIO.3,2X,5H EPM',EIO.3,IX,SH SM =,EIO.3, +3X,SH DF -,EIO.3,5X,5H ZQ ',EIO.3,7H SIZE =,EIO.3,//) FORMAT(2X, 6H CREF'.EIO.3,2X,AHEPSI,EIO.3,///) CONVERSION OF DT, ZQ, ALF, AND THE C-EXPANSION COEFFICIENTS TO DIMENSIONLESS QUANTITIES FOR USE IN THE MAIN PROGRAM CON - EPM/(2.0*96h85.0*ZP*ZM) DIM-CON/(A*GDE) DT - GDE*DT/CON ZQ - -CON*ZQ/(A*SM*EPM) ALF - CON*ALF/GDE DO 50 I-I,6 C(I) - (8.3IAE7)*TM*C(I)*CON/(A*A*GDE*2.0) CONTINUE C(2) - C(2)*(SM)**0-5 C(A) - C(h)*(SM)**O.5 6(5) 5 C(5)*(SM)**0.5 C(6) - c(6)*SM 6(7) . C(7)*SH C(8) - C(8)*SQRT(SM) C(9) - C(9)*SM C(IO) ' C(I0)*SM**2.0 DETERMINATION OF THE STEADY STATE ANALYTICAL SOLUTION. IF (BCOND .EQ. 0.0)CALL BCONDIT(II,ITOT,DOO,DA,SIZE) I2 I ITOT-II KDO'KDO*DIM KSO'KSO*DIM KDA-KDA*DIM KSA-KSA*DIM KDLFKDL*DIM KDR-KDR*DIM KSR-KSR*DIM KSL'KSL*DIM RETURN END nnnn nnnn I6 non nnnn 147 SUBROUTINE IANDW (Y,TIME,KOUNT,RCOND,XER,NPT,NAME,JNAME,M + .KFLAG,SIZE) IANDW CALCULATES THE POTENTIAL AS A FUNCTION OF X, THE TOTAL CHARGE AND MASS IN THE SYSTEM. AND WRITES THE OUTPUT FOR EACH TIME ROW SPECIFIED BY IJK. DIMENSION Y(3OO),XER(IOO),PSI(IOO).YER(IOO) DIMENSION E(IOO),YI(IOO),YJ(IOO) DIMENSION NPT(IO),JNAME(IO),NAME(IO) DIMENSION ZXARAY(IOO.3).ZARAY(IOO,3) DIMENSION PD(I50),PDT(I50) COMMON /ELEC/ EPS,ZI COMMON /PARAM/ KEY,R,Q,DT,F,TOL,RANGE.ISET COMMON/UNITS/SM,DF,GDE.EPM,A,ZP,ZM,CON,TM,DIM,CREF,ZREF COMMON /BOX/ IW,NPLT.IJK COMMON /PLTPTS/ POINT(IO) COMMON/CURRENT/II,I2 COMMON/DIFF/DI,Dz,DELTA COMMON/BNDRY/SL.SR,DL.DR INTEGER R,RR.Q.RSTOP,RSTART DOUBLE XER.XLO.XUP.AVINT,TMS,TCS,PSI,SIZE REAL II,I2 DATA I0/6)/.MAX/lQO/,ISCALE/l/,NF/3/ DATA ZXARAY/300*0.0/ RRI3*R RSTDP-R/z RSTART-RSTOP+I STRIPPING OUT THE MASS CHARGE.AND ELECTRIC FIELD VECTORS FOR INTEGRATION AND FORMING THE + AND - ION CONCENTRATIONS FOR PLOTS. DO I6 KII,R ZARAY(K,I)IZM*SM*(Y(3*K-l) - Y(3*K-2) + (OF/SM) - I.O)*I.0E-6 ZARAY(K,2)IZP*SM*(Y(3*K-I) + Y(3*K-2) + (OF/SM) + l.O)*I.0E-6 ZARAY(K,3) - Y(3*K)*I.0E-II*A*SM/CON YI(K) - Y(3*K-I) YJ(K) - Y(3*K-2) E(K) - Y(3*K) CONTINUE LOWER AND UPPER BOUNDS FOR INTEGRATION XLO I XER(I) XUP I XER(R) CALCULATING THE TOTAL CHARGE AND EXCESS MOLES OF IONS IN THE SYSTEM TCS-AVINT(XER,YI,R.XLO,XUP,IND) TMs- (ZP+ZM)*AVINT(XER.YI,R,XLO.XUP,IND) + + (ZP-ZM)*AVINT(XER,YJ,R,XLO,XUP,IND) no 26 23I nnnn 232 24 25 44A nnnnnnn 148 CALCULATING THE POTENTIAL AS A FUNCTION OF X DO 26 JII.R XUPIXER(J) PSI(J)I (I.OE-O7*A*SM/CON)*AVINT(XER,E,R,XLO,XUP,IND) CONTINUE RJ - PSI(R) ZJD - -D|M*(l2+ll)/(2.0*ZP*ZM*96485.0*SM) COMPUTING THE POTENTIAL LOSS DUE TO OHMIC DROP. PSOM - O.O PSOD - 0.0 no 231 J - 2,R SMM - Y(3*J-5) SPP - Y(3*J-2) DMM - Y(3*J-4) DPP - Y(3*J-I) S - O.5*(SMM + SPP) D = 0.5*(DMM + DPP) SETTING ARGUMENTS FOR THE FUNCTION SUBROUTINE. Tl~I ZIJ(S,D) SPAC - XER(J) - XER(J-l) PSOD - PSOD - (ZI+ZJD)*(I.0E-07*A*SM/CON)*SPAC/ZDE(S,D) PSOM - PSOD + PSOM CONTINUE PD(I+M/IJK) - PSI(R) - PSOM IF (M .EQ. 0 .OR. M .EQ. 1)GO TO 232 CHNG - (PD(I+M/|JK) - PD(I+(M-l)/|JK))/PD(I+M/IJK) ,CHNG - ABS(CHNG) IF (CHNG .LT. O.OOI)QIM CONVERTING EJ TO VOLTS, TIME TO SECONDS. TCS TO COULOMBS, TMS TO MOLES. AND ZI TO AMPS PER CM**2 RIME I TIME*CON/GDE PDT(I + M/IJK) I RIME RI I -A*SM*GDE*EPM*ZI/CON**2 TCS I TCS*2.0*ZP*ZM*96485.0*SM TMSISM*TMS WRITE (6I.24) RIME,KOUNT,TMS WRITE (6I.25) RJ,TCS.RI WRITE (6I.h44) RCOND FORMAT (IHI./.7X,IHX,IOX,3HSUM,9X,3HRHO.9X,2HC+,IOX,2HC-,IOX,IHE, IIOX, 3HPSI, 6X, 3H TI, EIO. 3, 2X, 7H KOUNTI, I2, 7X, 6H TMS I ,EIO. 3, /) FORMAT (86X, 4H EJI, EIO. 3, IX, 5H TCSI, EIO. 3, IX, 4H ZII, EIO. 3, /) FORMAT (86X, 8H RCOND I ,EIO. 3) AY IS THE SUM CONCENTRATION LESS THE BULK CONCENTRATION BY IS THE CHARGE DENSITY IN COULOMBS/CC CY IS THE ELECTRIC FIELD IN VOLTS/CM CP IS THE POSITIVE ION CONCENTRATION CM IS THE NEGATIVE ION CONCENTRATION 888 999 23 II IS 789 27 28 149 OD 888 KII.RST0P YER(K) - XER(K) + 0.5*SIZE CONTINUE DO 999 KIRSTART,R YER(K)I0.5*SIZE - XER(K) CONTINUE JR - RR — 2 DO 11 III.JR,3 AY - Y(l)*SM BY - Y(l+l)*SM*ZP*ZM*2.0*96485.0 CY - -Y(I+2)*l.OE-07*A*SM/CON CPIZM*((Y(I+I)+(DF/SM)) - (Y(I)+(I.o)))*SM CMIZP*((Y(I+I)+(DF/SM)) + (Y(I)+(I.O)))*SM WRITE (61,23) YER((I+2)/3).AY,BY,CP,CM,CY,PSI((I+2)/3) FORMAT (7E12.s) CONTINUE IF (M .EQ. Q) GO TO 15 RETURN MQ - I + Q/IJK NRITE (6I.28) NRITE (61,27) (PDT(I),PD(I).I-1,MQ) CPO - ZM*((Y(2)+(DF/SM)) - (Y(l)+(l.0)))*SM CMO - ZP*((Y(2)+(DF/SM)) + (Y(l)+(l.0)))*SM CPA - ZM*((Y(RR-l)+(DF/SM)) - (Y(RR-2)+(I.O)))*SM CMA - ZP*((Y(RR-l)+(DF/SM)) + (Y(RR-2)+(l.0)))*SM CPL - ZM*((DL+(DF/SM)) - (SL+(I.0)))*SM CML - ZP*((DL+(DF/SM)) + (SL+(I.0)))*SM CPR - ZM*((DR+(DF/SM)) - (SR+(I.0)))*SM CMR - ZP*((DR+(DF/SM)) + (SR+(l.0)))*SM ZIIL - (Dl*ZP*96485.0/DELTA)*(CPL-CPO) ZIZL - (02*2M*96A85.O/DELTA)*(CML-CMO) ZIIR - (DI*ZP*96A85.0/DELTA)*(CPA-CPR) ZIZR - (02*ZM*96A85.0/DELTA)*(CMA-CMR) NRITE (6I,789)ZIIL,ZI2L,ZIIR.ZI2R FORMAT (//,1x,5HIIL -,E12.5,3x,5H12L I,EI2.5,3X,5HIIR -,E12.5, + 3x,5HI2R -,E12.5.//) FORMAT (2X,2EI3.6) FORMAT (//.IX.46H POTENTIAL DIFFERENCE CORRECTED FOR OHMIC DROP./) RETURN END SUBROUTINE BCONDIT (II, ITOT, Doo, DA, SIZE) REAL 11,12,1TOT ‘COMMON /ELEC/ EPS.ZI COMMON/GREEK/C(2h) COMMON/UNITS/SM,DF,GDE.EPM,A,ZP.ZM,CON,TM,DIM,CREF,ZREF COMMON/THERM/QDELT,VI,V2 12 - ITOT-II n n on n 150 DETERMINATION OF STEADY STATE PARAMETERS ZJS - -DIM*(I2-II)/(2.0*ZP*ZM*96485.0*SM) ZJD - ~DIM*(l2+II)/(2.0*ZP*ZM*96485.0*SM) S - 0.0 D - 0.0 T1 - ZIJ(S,D) zz - (ZSSIS.D)*ZDD(S.D)) - (ZSD(S.D)*ZDS(S.D)) BDO - ((ZSS(S.D)*ZDE(S.D)) - (ZSEIS.D)*ZDS(S.D)))/ZZ BSO - ((ZDDIS.D)*ZSE(S.D)) - (ZSD(S.D)*ZDE(S.D)))/ZZ ..ADO - -((ZDS(S.D)*ZJ$) - (ZSS(S.D)*ZJD))/ZZ ASO - -((ZSD(S.D)*ZJD) - (ZDD(S.D)*ZJS))/ZZ AA - (-BDO/(A*A*EPS))**0.5 CALCULATION OF FIRST ORDER TERMS. BSD - ((ZDD(S.D)*ZSED(S.D)) - (ZSD(S.D)*ZDED(S.D)))/ZZ BSS - ((ZDD(S.D)*ZSEP(S.D)) - (ZSD(S.D)*ZDEP(S.D)))/ZZ BOD - ((ZSSIS.D)*ZDED(S.D)) - (ZDS(S.D)*ZSED(S.D)))/ZZ BDS - ((ZSS(S.D)*ZDEP(S.D)) - (ZDSIS.D)*ZSEP(S.D)))/ZZ ZZP - ZSS(S.D)*ZDDP(S.D) + ZSSP(S.D)*ZDD(S.D) + - ZSDIS.D)*ZDSP(S,D) - ZSDP(S,D)*ZDS(S,D) , zzn - ZSSIS.D)*ZDDD(S.D) + ZSSD(S.D)*ZDD(S.D) + - ZSDIS.D)*ZDSD(S.D) - ZSDD(S.D)*ZDS(S.D) TI - -((ZDSP(S.D)*ZJS) - (ZSSPIS.D)*ZJD)) T2 - -((ZDS(S.D)*ZJS) - (ZSSIS.D)*ZJD)) ' ADOP - (ZZ*TI - ZZP*T2)/(ZZ**2.0) TI - -((ZDSD(S.D)*ZJS) - (ZSSDIS.D)*ZJD)) ADOD - (22*Tl - ZZD*T2)/(ZZ**2.0) TI - -((ZSDP(S.D)*ZJD) - (ZDDP(S,D)*ZJS)) T2 - -((ZSD(S.D)*ZJD) - (ZDD(S.D)*ZJS)) ASOP - (ZZ*TI - ZZP*T2)/(ZZ**2.0) TI - -((ZSDD(S.D)*ZJD) - (ZDDD(S.D)*ZJS)) ASOD - (ZZ*Tl - zzo*T2)/(zz**2.o) CALCULATION OF THE THERMAL TERMS. YID+(DF/SM)-S-I.0 z- D + S + 1.0 + (OF/SM) IF (Y.GE.O.O) Y--1.OE-IO IF (z.LE.O.O) z-I.OE-IO SUM-((ZM*ZM)*Y - (ZP*ZP)*Z)/(ZM-ZP) ZPPIC(I)*(ZM*(Y)) + C(2)*(ZM*(Y))**I.5 Inn-C(3)*(ZP*(Z)) + C(h)*(ZP*(Z))**1-5 ZPMIC(5)*(SUM**I.5) + C(6)*(SUM**2.0) T1 - ZSSIS.D)*(ZP*ZPP+ZM*ZPM) + ZDS(S,D)*(ZP*ZPP-ZM*ZPM) T2 - ZSS(S,D)*(ZP*ZP*ZPP + 2.0*ZP*ZM*ZPM + ZM*ZM*ZMM) T3 - -ZDS(S.D)*(ZM*ZM*ZMM - ZP*ZP*ZPP) T5 - QDELT/(SIZE*VI) TA - Tl*T5*CON*BDO*l.0E+07/(A*96h85.0*h.18A*SM*(T2+T3)) ADO - ADO + TA THE FOLLONING ARE ELECTROLYTE CONCENTRATIONS AT THE NALLS, AND NOT RESERVOIR CONCENTRATIONS. n n n 1.78 III II2 II4 113 67 151 SAI - (((ASO*BDOI'(ADO*BSO))/BDOI*(SIZE/(2.0*A)) SA2-(BSO/BDo)*(DA-Doo)/2.o THE ZEROTH ORDER SUM CONCENTRATION AT THE WALLS. SRO - SAI + SA2 SLo--SRO ALPHA - ((ASO*BDO) - (ADO*BSO))/BDO GAMMA - SIZE/A SAI - (ALPHA*(GAMMA**2.0)/12.0)*(-ADO/BDO) SAI - SAI*(BSS-(BDS*BSO/BDO)) SA2 - ALPHA*GAMMA/2.0 SA2 - SA2*(BSS-(BDS*BSO/BDO)I*(DA/BDO) SA3 - 2.0*((DA/2.0)**2.0)/BDO SA3 - SA3*((BSD+(BSS*BSO/BDO)) - ((BSO/BDOI*(BDD+(BDS*BSO/BDO)))) THE FIRST ORDER CONTRIBUTION TO THE SUM CONCENTRATION AT THE NALLS. SR1 - SAI + SA2 + SA3 SL1 - SR1 SSR - SRO + SR1 SSL - SLO + SL1 COMPUTATION OF SECOND ORDER CONTRIBUTION.TO THE POTENTIAL DROP ACROSS THE ELECTROLYTE. T1 - (BDS*SIZE*ALPHA/(BDO*A))**2.0 T2 - -ADO*SIZE/(A*IZ.O*BDO) T3 - TI*T2 TA - 3.0*ALPHA*BDS/2.0 T5 - -BDS*BSO*ADO/(3.0*BDO) T6 - BDS*ALPHA*(DA-DOO)/(8.0*BDOI T7 - (SIZE/(A*BDO))**2.0 T8 - T6*T7*(TA-T5) T9 - T8 + T3 PHIZ - (l.0E-07)*A*A*SM*T9/CON NRITE (61,III)BDO,BSO,ADO,ASO NRITE (61.112)BSD.BSS,BDD,BDS NRITE (61,IIA)ADOP,ASOP,ADOD,ASOD THE VALUE OF PHI REPRESENTS THE ZEROTH ORDER CONTRIBUTION TO THE TOTAL POTENTIAL DIFFERENCE ACROSS THE CELL. PHI - -(I.OE-07*A*A*SM/(CON*BDOI)*(IDOO-DA) + ADO*SIZE/A) NRITE (61,113)AA,PHI NRITE(6I,67) 11,12 NRITE (61,555)SLO,SRO,SLI,SRI,SSL,SSR THE VALUE OF PHIz REPRESENTS THE SECOND ORDER CONTRIBUTION To THE POTENTIAL DROP ACROSS THE ELECTROLYTE. NRITE (61,A78)PHI2 FORMAT (2x,6HPH12 I,El0.3,//) FORMAT (IHI.Ix,5HBDO -,EIO.3,2x,5HBso -,EIO.3,2x,5HADO -,EIO.3, + 2X,5HASO I,EIO.3,//) FORMAT (2X,5HBSD I,EIO.3,2X,5HBSS I,EIO.3,2X,5HBDD I,EIO.3, + 2X,5HBDS I,EIO.3.//) FORMAT (2X,6HADOP I,EIO.3,2X,6HASOP I,EIO.3,2X,6HADOD '.EIO.3, + 2x,6HASOD -.E10.3.//) FORMAT (2X,4HAA I,EIO.3,3X,5HPHI I,EIO.3,5HVOLTS,//) FORMAT (2X,I5HIONIC CURRENT I,EI2.5,3X.20HELECTRONIC CURRENT I +.EIZ-5.//) nonnnnnn 555 80 20 9O 70 I00 30 40 152 FORMAT (2X,5HSLO '.EIO.3,3X,5HSRO I,EIO.3,3X,5HSLI I,EIO.3, + 3X,5HSRI I,EIO.3,3X,5HSSL I,EIO.3,3X,5HSSR I,EIO.3,//) RETURN END SUBROUTINE GRIDPR IDEL.XER.SIZE.DEBL) GRIDPR GENERATES A NONUNIFORM SPACE MESH WITH INCREASING DISTANCE BETWEEN GRID POINTS AS X INCREASES. IT IS DESIGNED FOR SIMULATIONS USING A POLARIZABLE ELECTRODE AT XIO AND A REVERSIBLE ELECTRODE AT XII. THE VALUE OF RANGE SHOULD BE SELECTED WITH CARE TO AVOID UNDERFLOW AND OVERFLOW FROM THE EXPONENTIAL FUNCTION DIMENSION DEL(IOO).XER(IOO) INTEGER R DOUBLE DEL,XER,SIZE.DEBL,GAUGE COMMON/UNITS/SM,DF,GDE.EPM,A,2P,2M,CON,TM,DIM,CREF,2REF COMMON /PARAM/ KEY,R,Q,DT,F,TOL,RANGE,ISET GAUGE - DEBL/5.O DO 10 I - I,ISET DEL(I) - GAUGE CONTINUE JSET - ISET + I L-MOD(R,2) IF (L.EQ.O)GO TO 70 JRI(R-I)/2 DO 20 JIJSET,JR DEL(J)IGAUGE*EXP((J-ISET)*RANGE) CONTINUE KSETIJR+2-L LSET-R-I DO 90 JIKSET.LSET DEL(J)IDEL(R-J) CONTINUE GO TO 100 JR-(R/2)-I DEL(R/ZIIGAUGE*EXP(((R/Z)-ISET)*RANGE) GO TO 80 JRIR-l SUM - 0.0 DO 30 M - I,JR SUM - SUM + DEL(M) CONTINUE DO A0 L - I,JR DEL(L) - DEL(L)*SIZE/(A*SUM) CONTINUE XER(I) - 0.0 no 50 M - 2,R nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnonnnnn 50 60 153 XER(M) - XER(M-l) + DEL(M-I) CONTINUE DO 60 NI2,JR XER(N) - A*XER(N) CONTINUE XER(R)-SIZE RETURN END SUBROUTINE COEMAT (R.DEL.DT.Y.ABD.M.MMM) SUBROUTINE COEMAT ASSEMBLES THE CORE MATRIX FOR USE IN THE NEWTON-RAPHSON ITERATION PROCEDURE. THE MATRIX IS ASSEMBLED DIRECTLY INTO BAND STORAGE FORM FOR USE IN THE LINPAC SUBROUTINE SGBCO. R IS THE NUMBER OF SPACE GRID POINTS. RR IS THE TOTAL NUMBER OF'DEPENDENT VARIABLES. DEL IS THE GRID SPACING VECTOR. DT IS THE TIME STEP. ABD(LDA,RR) IS THE MATRIX IN BAND STORAGE FORM. LDAIZ*ML+MU+I ML-NUMBER OF DIAGONALS ABOVE THE MAIN. MU-NUMBER OF DIAGONALS BELOW THE MAIN. Y IS THE SOLUTION VECTOR FROM THE PREVIOUS TIME ROW OR PREVIOUS ITERATION. COEMAT CALLS THE 2 COEFFICIENT FUNCTION SUBPROGRAMS. FOR THIS ROUTINE. THE CONVERSION FROM BAND FORM TO BAND STORAGE IS ABD(I-J+II.J)IA(I,J) THE BAND STORAGE MATRIX IS PASSED TO THE MAIN PROGRAM THROUGH THE COMMON BLOCK COE. THE RATE CONSTANTS ARE PASSED TO COEMAT THROUGH THE COMMON BLOCK RCNSTS. IF THERE WERE I5 SPACE GRID PTS. THE A AND ABD MATRICES WOULD LOOK LIKE THE FOLLOWING BLOCKS. A ABD c>c>c>c>c>c>c>c>c>uar-u>oo\10\ c>c>c:c>c>c>c>c>c>x-u>oo~JCh\n CDCDC)C)C>C>C)C>C>U>OD\JONU1£' c>c>c>c>c>c>u:x-uioo\IONUIr-u: c>c>c>c>c>c::-u>oo~so~v1r-uana CDC>CJCDCJCDUDOO\JO\U1£‘UJN)—‘ c>c>c:u:x-u>00\lonuwr-u:c>c>c> c>c>c>r-u:cn~4C>C) a:x-u>oo~10~uwz-u:c>c>c>c>c>c> r-uaoo~¢o~uwrruosac>c>c>c>c>c> Loc3c>c>c> oo~¢onvwrru¢c>c>c>c>c>c>c>c>c> \JONUTS'UPK!C>C)C>C>C>C>CDOHO o~uwr-uans—-c>c>c>c>c>c>c>c>c> a:r-u:oo~10~€>c>c>c>c>c>c><><3<3 c>x>uocn~4C>C>C)CDCDCJCD fill-U)OO\JO\U1S‘UJCDCJCDCDCDCDC) c:r-u:cm~4<>c>c>c>c> C)CDUJOO\JO\U1J‘UJh3—-C>C)C>C>C3 a!I'M)OD\JORU1JPMHCD(D(3C)C>C)C) u:r-u>oo~406U1srCo<3<3<3<3<3c>c>c>c>c> c>c>u3cnIJChunJruuhD-ac>c>c>c>c> c>c>c>anxioxuwsruac>c>c>c>c>c>c>. C)CDC)C)\JO\U1¥'UDK)C>C)C>C>CDC) c:c>c>c>c>oxuws-uan:-c>c>c>c>c> INTEGER R,RR,RSTOP nn nnnn n 100 154 REAL KDL,KSL,KDO,KSO,KDA,KSA,KDR.KSR DOUBLE DEL DIMENSION DEL(IOO),Y(300),ABD(I6,300) COMMON/UNITS/SM,DF,GDE.EPM.A.ZP.ZM.CON.TM.DIM,CREF,ZREF COMMON /ELEC/ EPS,ZI COMMON/RATES/KSL,KDL,KSO,KDO,KSA,KDA,KSR.KDR RR I 3*R RSTOP I RR - 5 INITIALIZING THE ABD MATRIX TO ZERO DO 100 lII.l6 00 100 JII.RR ABD(|.J) - 0.0 CONTINUE ASSEMBLING THE MATRIX ELEMENTS ARISING FROM THE BOUNDARY CONDITIONS AT XIO 02 . DEL(1)/2 O 022 - DEL(I)*DEL(I)/2.0 S - Y(I) D - Y(2) SP - Y(A) DP - Y(S) SETTING THE COEFFICIENT ARGUMENTS FOR THE FUNCTION ROUTINES. TI - ZIJ(S,D) T2 - ZIJUP(SP,DP) FI - 1.0/DT F2 - KSO/DZ F3 - ZSS(S.D)/022 TI - FI + F2 + F3 T2 - ZSSP(S.D)*Y(I)/DZZ F1 - ZSDP(S.D)/022 F2 - -ZSEP(S.D)/EPS T3 - (F1 + F2)*Y(2) TA - ZSEP(S.D)*Y(3)/(2.0*DZ) T5 - -ZSSP(S.D)*Y(A)/Dzz T6 - -ZSDP(S.D)*Y(5)/DZZ ABD(II.I) I Tl + T2 + T3 + Th + T5 + T6 FI - KOO/DZ F2 - ZDS(S.D)/022 T1 - F1 + F2 T2 - ZDSP(S.D)*Y(I)/022 FI - ZDDP(s.D)/022 F2 - -ZDEP(S.D)/EPS T3 - (FI + F2)*Y(2) TA ZDEP(S,D)*Y(3)/(2.0*02) T5 - -ZDSP(S.D)*Y(h)/022 T6 - -ZDDP(S,D)*Y(5)/022 nonnn 155 ABD(12,1) - Tl + T2 + T3 + TA + T5 + T6 ABD(I3,I) I KDO/EPS TI - KDO/DZ' T2 - ZSD(S,D)/022 T3 - -ZSE(S.D)/EPS Th - ZSSD(S.D)*Y(I)/DZZ T5 . ZSDDIS.D)*Y(2)/022 T6 - -ZSSD(S.D)*Y(A)/022 T7 - -ZSDD(S,D)*Y(5)/DZZ T8 - -ZSED(S.D)*Y(2)/EPS T9 - ZSED(S,D)*Y(3)/(2.0*DZ) . ABD(IO,2) . T1+T2+T3+TA+T5+T6+T7+T8+T9 TI - 1.0/DT T2 - KSO/DZ T3 - ZDD(S.D)/022 Th - -ZDE(S.D)/EPS T5 - ZDSD(S.D)*Y(I)/022 T6 - ZDDDIS.D)*Y(2)/022 T7 - -ZDSD(S.D)*Y(4)/022 T8 - -ZDDD(S.D)*Y(5)/DZZ T9 - -ZDED(S.D)*Y(2)/EPS T10 - ZDED(S,D)*Y(3)/(2.0*DZ) ABD(II,2) - TI+T2+T3+TA+T5+T6+T7+T8+T9+T10 ABD(12.2) - KSO/EPS TI - ZSEIS.D)/(2.0*DZ) T2 - ZSEUP(SP,DP)/(2.0*DZ) ABD(9.3) - TI + T2 T1 - ZDE(S.D)/(2.0*DZ) T2 - ZDEUP(SP,DP)/(2.0*DZ) ABD(10,3) - TI + T2 ABD(II,3) - 1.0/DT ABD(8,A) - -ZSS(S,D)/022 + (0.5*ZSEPUP(SP.DP)/DZ)*Y(3) ABD(9,A) - -ZDS(S,0)/022 + (O.5*ZDEPUP(SP,DP)/DZ)*Y(3) ABD(10,A) - 0.0 ABD(7.5) - -ZSD(S,D)/022 + (0.5*ZSEDUP(SP,DP)/DZ)*Y(3) ABD(8,5) - -ZDD(S,D)/022 + (0.5*ZDEDUP(SP,DP)/DZ)*Y(3) ABD(9.5) - 0.0 ABD(6.6) - 0.0 ABD(7.6) - 0.0 ABD(8.6) - 0.0 ASSEMBLING THE INTERIOR MATRIX ELEMENTS IT IS CONVENIENT TO LOOP IN MULTIPLES OF 3 BECAUSE THERE ARE 3 EQUATIONS AT EVERY GRID POINT DO 101 K-A,RST0P,3 KM - (K-1)/3 KK - (K+2)/3 SF - (Y(K-3) + Y(K))/2.0 SMM - Y(K-3) SPP . Y(K+3) S - Y(K) SP - (Y(K+3) + Y(K))/2.0 n 156 on - (Y(K-z) + Y(K+I))/2.0 DMM - Y(K—z) D - Y(K+l) DP - (Y(K+A) + Y(K+l))/2.0 OPP - Y(K+A) AD - DEL(KM) + DEL(KK) DENM - DEL(KM)*AD/2.0 DENK - DEL(KK)*AD/2.0 RAT - DEL(KK)/DEL(KM) . AMUL I DEL(KK)*DEL(KM) SUB - DEL(KK) - DEL(KM) SETTING THE COEFFICIENT ARGUMENTS FOR THE FUNCTION ROUTINES. TI I ZIJ(S,D) T2 I ZIJUPISP,DP) T3 - ZIJDN(SF,DM) TA - PPZIJ(SPP,DPP) T5 - DDZIJ(SMM,DMM) T1 - —ZSSDN(SF.DM)/DENM T2 - -ZSSPDN(SF.DM)*Y(K-3I/DENM T3 - -ZSDPDN(SF,DM)*Y(K-2)/DENM TA . ZSSPDN(SF,DM)*Y(K)/DENM T5 - ZSDPDN(SF.DM)*Y(K+I)/DENM T6 - -DDZSEP(SMM,DMM)*RAT*Y(K+2)/AD ABD(IA,K-3) a T1 + T2 + T3 + TA + T5 + T6 T1 - —ZDSDN(SF,DM)/DENM T2 - -ZDSPDN(SF.DM)*Y(K'3)/DENM T3 - -ZDDPDN(SF.DM)*Y(K-2)/DENM TA - ZDSPDN(SF.DM)*Y(K)/DENM T5 - ZDDPDN(SF,DM)*Y(K+l)/DENM T6 - -DDZDEP(SMM,DMM)*RAT*Y(K+2)/AD ABD(15,K-3) - TI + T2 + T3 + TA + T5 + T6 ABD(I6,K-3) - -ZDS(S.D)*RAT/(EPS*AD) T1 - -ZSDDN(SF.DM)/DENM T2 - -ZSSDDN(SF.DM)*Y(K-3)/DENM T3 - -ZSDDDN(SF,DM)*Y(K-2)/DENM TA - ZSSDDN(SF,DM)*Y(K)/DENM T5 - ZSDDDN(SF,DM)*Y(K+l)/DENM T6 - -DDZSED(SMM.DMM)*RAT*Y(K+2)/AD ABD(I3,K-2) - T1 + T2 + T3 + TA + T5 + T6 T1 - -ZDDDN(SF,DM)/DENM T2 - -ZDSDDN(SF.DM)*Y(K-3)/DENM T3 - -ZDDDDN(SF.DM)*Y(K-ZI/DENM TA - ZDSDDN(SF,DM)*Y(K)/DENM T5 . ZDDDDN(SF,DM)*Y(K+l)/DENM T6 - -DDZDEDISMM.DMM)*Y(K+2)*RAT/AD ABD(IA,K-2) T1 + T2 + T3 + TA + T5 + T6 ABD(I5,K-2) -ZDD(S,D)*RAT/(EPS*AD) ABD(12,K-1) . 0.0 ABD(I3,K-l) - 0.0 ABD(Ith'I) . 0.0 157 TI I I.O/DT T2 I ZSSUP(SP,DP)/DENK T3 I ZSSDNISF.DM)/DENM Th I -ZSEP(S,D)*Y(K+I)/EPS T5 I ZSEP(S,D)*SUB*Y(K+2)/AMUL ABD(II,K) . T1 + T2 + T3 + TA + T5 T1 - ZDSUP(SP.DP)/DENK T2 - ZDSDN(SF,DM)/DENM T3 - -ZDEP(S.D)*Y(K+l)/EPS TA - ZDEP(S.D)*Y(K+2)*SUB/AMUL ABD(12.K) I TI + T2 + T3 + TA TI - -ZDSP(S,D)*Y(K-3)*RAT/(EPS*AD) T2 - -ZDDP(S.D)*Y(K-2)*RAT/(EPS*AD) T3 - ZDSP(S.D)*Y(K)*SUB/(EPS*AMUL) TA - ZDDPIS.D)*Y(K+l)*SUB/(EPS*AMUL) T5 - -ZDEP(S.D)*Y(K+2)/EPS T6 - ZDSPIS.D)*Y(K+3)/(RAT*EPS*AD) T7 - ZDDPIS.D)*Y(K+h)/(RAT*EPS*AD) T8 - ZDSIS.D)*SUB/(EPS*AMUL) ABD(I3,K) - TI + T2 + T3 + TA + T5 + T6 + T7 + T8 T1 - ZSDUPISP,DP)/DENK T2 - ZSDDN(SF,DM)/DENM T3 - -ZSE(S.D)/EPS TA - -ZSED(S.D)*Y(K+I) T5 - ZSED(S.D)*SUB*Y(K+2)/AMUL ABD(IO.K+I) - TI + T2 + T3 + Th + T5 TI I I.O/DT T2 I ZDDUP(SP,DP)/DENK T3 ' ZDDDN(SF.DM)/DENM TA I ’ZDEIS.D)/EPS T5 I -ZDED(S,D)*Y(K+I) T6 - ZDEDIS.D)*SUB*Y(K+2)/AMUL ABD(II,K+I) - TI + T2 + T3 + TA + T5 + T6 Tl - -ZDSD(S.D)*Y(K-3)*RAT/(EPS*AD) T2 - ~ZDDD(5.0)*Y(K-2)*RAT/(EPS*AD) T3 - ZDSD(S,D)*Y(K)*SUB/(EPS*AMUL) TA - ZDDDIS.D)*Y(K+I)*SUB/(EPS*AMUL) T5 - -ZDED(S.D)*Y(K+2)/EPS T6 - ZDSD(S.D)*Y(K+3)/(RAT*EPS*AD) T7 - ZDDD(S.D)*Y(K+A)/(RAT*EPS*AD) T8 - ZDD(S.D)*SUB/(EPS*AMUL) ABD(12,K+I) - T1 + T2 + T3 + TA + T5 + T6 + T7 + T8 T1 - -DDZSE(SMM,DMM)*RAT/AD T2 - PPZSE(SPP,DPP)/(RAT*AD) T3 - ZSEIS.D)*SUB/AMUL ABD(9,K+2) - TI + T2 + T3 T1 - -DDZDE(SMM.DMM)*RAT/AD T2 - PPZDE(SPP,DPP)/(RAT*AD) T3 - ZDE(S.D)*SUB/AMUL ABD(IO,K+2) - Tl + T2 + T3 T1 - I.O/DT T2 - -ZDE(S.D)/EPS ABD(II.K+2) - TI + T2 nnnn no IOI 158 TI - -ZSSUPISP.DP)/DENK T2 - ZSSPUP(SP,DP)*Y(K)/DENK T3 - ZSDPUP(SP,DP)*Y(K+l)/DENK TA - -ZSSPUP(SP,DP)*Y(K+3)/DENK T5 - -ZSDPUP(SP,DP)*Y(K+A)/DENK T6 - PPZSEP(SPP,DPP)*Y(K+2)/(AD*RAT) ABD(8,K+3) - TI + T2 + T3 + TA + T5 + T6 T1 - -ZDSUP(SP.DP)/DENK T2 - ZDSPUP(SP,DP)*Y(K)/DENK T3 - ZDDPUP(SP.DP)*Y(K+l)/DENK TA - -ZDSPUP(SP,DP)*Y(K+3)/DENK T5 - -ZDDPUP(SP,DP)*Y(K+A)/DENK T6 - PPZOEP(SPP,DPP)*Y(K+2)/(AD*RAT) ABD(9,K+3) - Tl + T2 + T3 + TA + T5 + T6 ABD(IO,K+3) - ZDS(S.DI/(EPS*RAT*AD) TI - -ZSDUP(SP.DP)/DENK T2 - ZSSDUP(SP,DP)*Y(K)/DENK T3 - ZSDDUP(SP,DP)*Y(K+I)/DENK TA - -ZSSDUP(SP,DP)*Y(K+3)/DENK T5 - -ZSDDUP(SP,DP)*Y(K+4)/DENK T6 I PPZSED(SPP,DPP)*Y(K+2)/(AD*RAT) ABD(7,K+A) I TI + T2 + T3 + Th + T5 + T6 TI I -ZDDUP(SP,DP)/DENK T2 I ZDSDUP(SP,DP)*Y(K)/DENK T3 - ZDDDUPISP.DP)*Y(K+l)/DENK TA - -ZDSDUP(SP,DP)*Y(K+3)/DENK T5 . -ZDDDUP(SP,DP)*Y(K+4)/DENK T6 - PPZDED(SPP,DPP)*Y(K+2)/(AD*RAT) ABD(8,K+A) - T1 + T2 + T3 + TA + T5 + T6 ABD(9,K+A) ZDD(S.D)/(EPS*RAT*AD) ABD(7,K+5) - 0.0 ABD(8,K+5) I 0.0 CONTINUE ASSEMBLING THE MATRIX ELEMENTS ARISING FROM THE BOUNDARY CONDITIONS AT XIA 02 - DEL(R-I)/2.0 022 - DEL(R-l)*DEL(R-l)/2.0 S - Y(RR-z) SF - Y(RR-S) D - Y(RR-I) DM - Y(RR-A) SETTING THE COEFFICIENT ARGUMENTS FOR THE FUNCTION ROUTINES. TI - ZIJ(S.D) T2 - ZIJDN(SF.DM) ABD(IA,RR-5) I -ZSS(S,D)/DZZ - (0.5*ZSEPDN(SF.DM)/DZ)*Y(RR) ABD(I5,RR-5) I -ZDS(S,D)/022 - (O.5*ZDEPDN(SF.DM)/DZ)*YIRR) ABD(I6,RR-5) I 0.0 ABD(I3.RR-A) ABD(IA.RR-A) ABD(15,RR-A) ABD(12.RR-3) ABD(I3.RR-3) ABD(IA,RR-3) F1 - F2 - F3 T1 T2 T3 TA F1 F2 T5 T6 159 0.0 GOO GOO I.O/DT KSA/DZ zss(s.D)/022 FI + F2 + F3 -ZSSP(S,D)*Y(RR-5)/022 -ZSDP(S,D)*Y(RR-A)/022 ZSSPIS.D)*Y(RR-2)/022 ZSDP(S,D)/022 -ZSEP(S,D)/EPS (F1 + F2)*Y(RR-1) -ZSEP(S.D)*Y(RR)/(2.0*DZ) ABD(II,RR-2) I T1 + T2 + T3 + TA + T5 + T6 FI I F2 I TI I T2 T3 TA FI F2 T5 T6I KDA/02 ZDS(S.D)/022 F1 + F2 -ZDEP(S,D)*Y(RRI/(2.0*02) -ZDDP(S,D)*Y(RR-A)/022 ZDSP(S.D)*Y(RR-2)/022 ZDDP(S,D)/022 -z0EP(S.D)/EPS (FI + F2)*Y(RR- 1) -ZDSP(S, D)*Y(RR- 5)/022 ABD(IZ. RR- 2) I T) + T2 + T3 + TA + T5 + T6 ABD(IL RR- 2) I -KDA/EPS T1 T2 T3 TA T5 T6 T7 T8 T9- IKDA/DZ ZSD(S, D)/022 -ZSE(S.D)/EPS -ZSSD(S,D)*Y(RR-5)/022 -ZSDD(S,D)*Y(RR-A)/022 ZSSD(S.D)*Y(RR-2)/022 ZSDD(S.D)*Y(RR-I)/022 -ZSED(S. D)*Y(RR- I)/EPS -ZSED(S, D)*Y(RR)/(2. 0*02) ABD(IO, RR- 1) I TI+T2+T3+TA+T5+T6+T7+T8+T9 TI - T2 - T3 TA T5 T6 T7 - T8 - T9 - I. O/DT KSA/DZ 200(S.D)/022 -z0E(S.D)/EPS -ZDSD(S.D)*Y(RR-5)/022 -ZDDD(S.D)*Y(RR-A)/022 ZDSD(S.D)*Y(RR-2)/022 ZDDD(S,D)*Y(RR-l)/022 -ZDED(S.D)*Y(RR-I)/EPS TIO I -ZDED(S.D)*Y(RR)/(2.0*DZ) ABD(11,RR-1) . T1+T2+T3+TA+T5+T6+T7+T8+T9+T10 ABD(12,RR-I) I -KSA/EPS ABD(9,RR) - -ZSE(S,D)/(2.0*02) - ZSEDN(SF.DM)/(2.0*DZ) -ZSD(S,D)/022 - (O.5*ZSEDDN(SF,DM)/DZ)*Y(RR) IZDD(S,D)/DZZ - (0.5*ZDEDDN(SF,DM)/DZ)*Y(RR) nnnnnnnnnnnnnnnn an n I02 I03 160 ABD(10,RR) - -ZDE(S.D)/(2.0*DZ) - ZDEDN(SF.DM)/(2.0*DZ) ABD(11,RR) - I.O/DT IF (M .NE. MMM) GO TO 103 OD 102 IM - 1.16 00 102 JM - 1.RR NRITE (61,1) ABD(IM,JM) FORMAT (2x,EIA.7) CONTINUE CONTINUE RETURN END SUBROUTINE RHSVN (R,DEL,DT,Y,B) WHEN IFLAGIO (THE FIRST TIME RHSV IS CALLED AT A GIVEN TIME ROW) SUBROUTINE RHSV COMPUTES THE RHS VECTOR FOR THE MATRIX EQUATION AT THE N+I TH TIME ROW FROM THE SOLUTION TO THE MATRIX EQUATION AT THE N TH TIME ROW. WHEN IFLAGII. RHSV COMPUTES THE CURRENT VALUE OF THE AUXILLARY VECTOR. IN THE MAIN PROGRAM THIS VECTOR WILL BE ADDED TO THE RHS VECTOR TO FORM THE CURRENT RHS VECTOR FOR A GIVEN ITERATION. R IS THE NUMBER OF SPACE PTS. RR IS THE TOTAL NUMBER OF DEPENDENT VARIABLES. DT IS THE TIME STEP SIZE. DEL IS THE NON-UNIFORM GRID SPACING VECTOR. Y(RR) IS THE SOLUTION VECTOR FROM THE N TH TIME ROW IF IFLAGIO AND THE SOLUTION TO THE PREVIOUS ITERATION STEP IF IFLAGII. BIRR) IS THE RIGHT HAND SIDE VECTOR IF IFLAGIO AND THE AUXILLARY VECTOR IF IFLAGII. THE ONLY DIFFERENCES BETWEEN THE RHS VECTOR AND THE AUX VECTOR ARE IN THE TERMS INVOLVING DT AND IN THE BOUNDARY TERMS. INTEGER R.RR.RSTOP REAL KDL,KSL,KDO,KSO,KDA,KSA.KDR,KSR DOUBLE DEL DIMENSION DEL(IOO),Y(300),B(300) COMMON/UNITS/SM,DF,GDE.EPM,A,ZP,ZM,CON,TM,DIM,CREF,ZREF COMMON /ELEC/ EPS.ZI COMMON/RATES/KSL,KDL,KSO,KDO,KSA,KDA,KSR,KDR RR I 3*R RSTOP I RR - 5 VECTOR ELEMENTS ARISING FROM THE BOUNDARY CONDITIONS AT XIO 02 - DEL(I)/2.0 022 - DEL(I)*DEL(l)/2.o s - Y(I) SP - Y(A) D - Y(2) DP - Y(5) SETTING THE COEFFICIENT ARGUMENTS FOR THE FUNCTION ROUTINES. nnnnn n 161 TI I ZIJ(S,D) T2 - ZIJUP(SP.DP) F1 -I.O/DT F3 -KSO/DZ TI (F1 + F2 + F3)*Y(I) T2 (-ZSD(S.D)/022 + ZSE(S.D)/EPS - KDO/DZ)*Y(2) T3 -(ZSE(S.D)/(2.0*DZ) + ZSEUP(SP.DP)/(2.0*DZ))*Y(3) TA (ZSS(S.D)/022)*Y(A) T5 (ZSD(S.D)/022)*Y(5) B(1 - TI + T2 + T3 + TA + T5 IIIIINIIVIIIIIIII T1 (-KDO/Dz - ZDS(S.D)/DZZ)*Y(I) F1 -I.0/0T F2 -ZOD(S.D)/022 F3 ZDE(S.D)/EPS FA -KSO/DZ T2 (F1 + F2 + F3 + FA)*Y(2) T3 ~(ZDEIS.D)/(2.0*02) + ZDEUP(SP.DP)/(2.0*DZ))*Y(3) TA (ZDSIS.D)/022)*Y(A) T5 I (ZDD(S,D)/022)*Y(5) 8(2) - T1 + T2 + T3 + TA + T5 T1 - (-I.O/DT)*Y(3) T2 - -KSO*Y(2)/EPS T3 I -KDO*Y(l)/EPS 8(3) - TI + T2 + T3 VECTOR ELEMENTS FROM THE INTERIOR SPACE POINTS THE DO-LOOP IS IN MULTIPLES OF 3 SINCE THERE ARE 3 EQUATIONS AT EACH SPACE POINT 00 101 K-A,RST0P.3 KM I (K-II/3 KK . (K+2)/3 SF . (Y(K-3) + Y(K))/2.0 SMM - Y(K-3) SPP - Y(K+3) S - Y(K) SP - (Y(K+3) + Y(K))/2.0 DM - (Y(K-2) + Y(K+I))/2.0 DMM - Y(K-z) D . Y(K+1) DP.I (Y(K+A) + Y(K+l))/2.0 DPP - Y(K+A) AD - DEL(KM) + DEL(KK) DENM - DEL(KM)*AD/2.0 DENK - DEL(KK)*AD/2.0 RAT - DEL(KK)/DEL(KM) AMUL - DEL(KK)*DEL(KM) SUB - DEL(KK) - DEL(KM) SETTING THE COEFFICIENT ARGUMENTS FOR THE FUNCTION ROUTINES. n n 162 TI - ZIJ(S,D) T2 - ZIJUP(SP,DP) T3 - ZIJDN(SF,DM) TA - PPZIJ(SPP,DPP) T5 - ODZ|J(SMM,DMM) T1 - ZSSDN(SF,DM)*Y(K-3)/DENM T2 - ZSDONISF.DM)*Y(K-2)/DENM FI - -DOZSE(SMM.DMM)*RAT/AD F2 - PPZSE(SPP,DPP)/(RAT*AD) F3 - ZSE(S.D)*SUB/AMUL T3 - -(FI + F2 + F3)*Y(K+2) F1 - -1.0/DT F2 - -ZSSUP(SP,OP)/DENK F3 - -ZSSON(SF,DM)/OENM TA - (F1 + F2 + F3)*Y(K) F1 - -ZSOUP(SP.DP)/DENK F2 - -zs00N(SF.DM)/DENM F3 - ZSE(S.D)/EPS T5 . (F1 + F2 + F3)*Y(K+I) T6 - ZSSUP(SP,DP)*Y(K+3)/DENK T7 - ZSDUP(SP,DP)*Y(K+A)/DENK VECTOR ELEMENTS FROM S EQUATION 8(K) . T1 + T2 + T3 + TA + T5 + T6 + T7 Tl - ZDSON(SF,DM)*Y(K-3)/DENM T2 - -ZDSDN(SF.DM)*Y(K)/DENM T3 - ZDDON(SF,DM)*Y(K-2)/DENM TA - -ZDOON(SF.DM)*Y(K+l)/DENM T5 - ZDDUP(SP,DP)*Y(K+A)/DENK T6 - -ZDOUP(SP,DP)*Y(K+l)/DENK FI - -1.0/DT F2 I ZDE(S.D)/EPS T7 . (F1 + F2)*Y(K+I) T8 - ZOSUP(SP,DP)*Y(K+3)/DENK T9 - -ZDSUP(SP.DP)*Y(K)/DENK F1 - -DDZDE(SMM.DMM)*RAT/AO F2 - PPZDEISPP,DPP)/(RAT*AO) F3 I ZDE(S,D)*SUB/AMUL TIO I -(FI + F2 + F3)*Y(K+2) VECTOR ELEMENTS FROM D EQUATION B(K+I) - T1 + T2 + T3 + TA + T5 + T6 + T7 + T8 + T9 + T10 T1 - z05(s, D)*Y(K- 3)*RAT/(EPS*AO) T2 - z00(s, D)*Y(K- 2)*RAT/(EPS*AD) T3 - -ZDS(5.0)*YIK)*SUB/(EPS*AMUL) TA - -ZDO(5.0)*Y(K+l)*SUB/(EPS*AMUL) T5 - (-I.0/DT + ZDE(S,D)/EPS)*Y(K+2) T6 - 'ZDS(S.D)*Y(K+3)/(EPS*RAT*AD) T7 -ZDD($.D)*Y(K+A)/(EPS*RAT*AD) n on no IOI 163 VECTOR ELEMENTS FROM E EQUATION B(K+2) - TI + T2 + T3 + TA + T5 + T6 + T7 CONTINUE ~ VECTOR ELEMENTS ARISING FROM THE BOUNDARY CONDITIONS AT XIA 02 - DEL(R-I)/2.0 022 - DEL(R-I)*DEL(R-I)/2.0 s - Y(RR-z) SF - Y(RR-5) D - Y(RR-I) DM - Y(RR-A) SETTING THE COEFFICIENT ARGUMENTS FOR THE FUNCTION ROUTINES. TI I ZIJ(S,D) T2 I ZIJDN(SF.DM) TI - ZSS(S.D)*Y(RR-5)/DZZ T2 - ZSD(S.D)*Y(RR-A)/022 F1 - -I.0/DT F2 - -KSA/DZ F3 - -zss(s.D)/022 T3 - (F1 + F2 + F3)*Y(RR-2) FI - -KDA/D2 F2 - -ZSD(S.D)/022 F3 - ZSE(S.D)/EPS TA - (F1 + F2 + F3)*Y(RR-I) T5 - (ZSE(S,D)/(2.0*DZ) + ZSEDN(SF,DM)/(2.0*DZ))*Y(RR) B(RR—z) - T1 + T2 + T3 + TA + T5 T1 - ZDS(S.D)*Y(RR-5)/DZZ T2 - ZDD(S,D)*Y(RR-A)/022 T3 I (-KDA/D2 -ZDS(S,D)/DZZ)*Y(RR-2) FI I II.O/DT F2 I -KSA/DZ FA - ZDE(S.D)/EPS Th I (Fl + F2 + F3 + FA)*Y(RR-I) T5 I (ZDE(S,D)/(2.0*DZ) + ZDEDN(SF,DM)/(2.0*DZ))*Y(RR) B(RR-l) I TI + T2 + T3 + TA + T5 Tl - (-I.O/DT)*Y(RR) T2 I KSA*Y(RR-I)/EPS T3 I KDA*Y(RR-2)/EPS B(RR) I T] + T2 + T3 RETURN END nnnn I00 I09 164 SUBROUTINE RHSVO (R.DEL.DT.Y.B.M.MMM) SUBROUTINE RHSVO GENERATES THE PORTION OF THE B VECTOR THAT DOES NOT CHANGE WITH EACH ITERATION. INTEGER R,RR,RSTOP REAL KDL,KSL,KDO,KSO,KDA,KSA,KDR,KSR REAL II,I2 . DOUBLE DEL DIMENSION DEL(IOO).Y(300),B(300) COMMON /ELEC/ EPS.ZI COMMON/UNITS/SM,DF,GDE.EPM.A,ZP,ZM,CON,TM,DIM,CREF,ZREF COMMON/RATES/KSL,KDL,KSO,KDO.KSA.KDA,KSR,KDR COMMON/CURRENT/II,I2 COMMON/BNDRY/SL.SR.DL.DR . RR I 3*R RSTOP-RR-3 ZJD - -DIM*(Il+|2)/(2.0*ZP*ZM*96A85.0*SM) TlIY(I)/DT T2 - 2.0*(KSL*(SL+I.O) + KDL*(DL+(DF/SM)))/DEL(I) T3 - -2.0*(KSO + KDO*(DF/SM))/DEL(I) 8(1)-T1+T2+T3 TI-Y(2)/DT T2 - 2.0*(KSL*(DL+(DF/SM)) + KDL*(SL+I.O))/DEL(I) T3 - -2.0*(KSO*(DF/SM) + K00)/DEL(I) 8(2)ITI+T2+T3 TI-Y(3)/DT + z1/EPS T2 - ZJD/EPS T3 - (KSL*(DL+(DF/SM)) + KDL*(SL+I.0))/EPS TA - -(KSO*(DF/SM))/EPS - (KDO/EPs) 8(3)-TI+T2+T3+TA DO 100 I-A,RST0P B(I) - Y(I)/DT LIMOD(I,3) IF (L .EQ. 0) B(l)IB(I)+Zl/EPS + ZJD/EPS CONTINUE TI-Y(RR-2)/DT T2 - 2.0*(KSR*(SR+I.O) + KDR*(DR+(DF/SM)))/DEL(R-l) T3 - -2.0*(KSA + KDA*(DF/SM))/DEL(R-l) B(RR-2)-TI+T2+T3 TI-Y(RR-1)/DT T2 - 2.0*(KSR*(DR+(DF/SM)) + KDR*(SR+I.O))/DEL(R-l) T3 - -2.0*(KSA*(DF/SM) + KDA)/DEL(R-l) B(RR-I)-TI+T2+T3 T1-Y(RR)/OT + ZI/EPS T2 - ZJD/EPS T3 - -(KSR*(DR+(DF/SM)) + KDR*(SR+I.0))/EPS TA - (KSA*(DF/SM) + KDAI/EPS B(RR)-TI+T2+T3+TA IF (M .NE. MMM)GO TO 110 00 109 I-I,RR NRITE (61,1) B(I),I CONTINUE nonnnnnnnnnnnn an an no on n IIO I 100 IOI 165 CONTINUE FORMAT (2X,EIA.6,I5) RETURN END SUBROUTINE CMPR (Y.YY,TOL,R,KFLAG) SUBROUTINE CMPR COMPARES THE CURRENT ITERATION VECTOR, Y, TO THE CURRENT SOLUTION VECTOR, YY+Y. (YY+Y IS RENAMED YY DESTROYING THE OLD SOLUTION VECTOR) THE CURRENT ITERATION VECTOR IS TESTED TO DETERMINE IF IT IS SMALL ENOUGH FOR THE PROGRAM TO PROCEED TO THE NEXT TIME STEP. SINCE SOME OF THE SOLUTION VECTOR COMPONENTS ARE MUCH SMALLER AND HENCE LESS ACCURATE THEN OTHERS, THEY ARE NOT INCLUDED IN THE TESTING. EACH OF THE 3 COMPONENTS MAKING UP THE SOLUTION VECTOR, S,D, AND E ARE CHECKED SEPARATELY. KFLAGIO MEANS THE PRESENT ITERATION VECTOR IS SMALL ENOUGH FOR THE PROGRAM TO PROCEED TO THE NEXT TIME STEP. KFLAGII MEANS THE PROGRAM WILL ITERATE AGAIN WITH THE DERIVATIVES EVALUATED AT YY+Y. TOL IS THE RELATIVE CHANGE TOLERANCE ALLOWED. DIMENSION Y(300),YY(300) INTEGER R.RR KFLAGIO RRI3*R ‘FORMING THE NEN SOLUTION VECTOR DO 100 III,RR YY(I)IYY(I)+Y(I) CONTINUE M-I STRIPS OUT THE S COMPONENTS, M=2 THE D COMPONENTS, M=3 THE E COMPONENTS. DO 103 MIl,3 YMAx-ABS(YY(M)) FINDING THE MAXIMUM VALUE OF 5 (IF M-I),D (IF M-2), 0R E (IF M-3) DO 101 I-M,RR,3 YMAx-AMAx1(ABS(YY(I)),YMAx) CONTINUE FINDING THE MINIMUM VALUE THAT NILL BE CHECKED YMIN I YMAX/I.0E+09 D0 I02 JIM.RR.3 TESTING TO DETERMINE IF A VALUE SHOULD BE CHECKED n nonnn nnnnnnn I02 I03 IOA I0 166 IF (ABS(YYIJ)) .LT. YMIN) GO TO I02 CHECKING TO DETERMINE IF THE CHANGE IN THE SOLUTION VECTOR IS SMALL ENOUGH AT A GIVEN J THE FOLLONING 2 STATEMENTS PREVENT DIVISION BY ZERO AYIABS(YY(J)) IF (AY .LT. I.OE-99) GO TO 102 RATIOIY(J)/YY(J). TEST-ABS(RATIO) IF (TEST .GT. TOL) GO TO 10A CONTINUE CONTINUE RETURN KFLAG-I RETURN END FUNCTION AVINT (X,Y,N,XLO,XUP,IND) AVINT IS USED TO INTEGRATE THE NONUNIFORM SPACE MESH FOR COMPUTATION OF THE POTENTIAL DIFFERENCE, TOTAL CHARGE IN THE SYSTEM. AND TOTAL MOLES IN THE SYSTEM IF IND IS RETURNED AS -I XUP IS LESS THAT XLO AND THE COMPUTATION WILL NOT PROCEED DIMENSION X(N),Y(N) DOUBLE X,XLO,XUP,SYL,A,B,C,CA,CB,CC,TERMI,TERM2,TERM3 DOUBLE XI,X2,X3,SUM.DIFF INDIO IF (N .LT. 3) RETURN DO I0 II2.N IF (X(I) .LE. X(I-I)) RETURN CONTINUE SUM-0.0 IF (XLO .LE. XUP) GO TO 5 SYLIXUP XUPIXLO XLOISYL INDI-I GO TO 6 INDII SYLIXLO IBII JIN DO I III,N IF (X(I) .GE. XLO) GO TO 7 IBIIB+I con IA IS 167 IB-MAX0(2,IB) IB-MIN0(IB,N-I) DO 2 III,N IF (XUP .GE. X(J)) GO TO 8 JIJ-l JIMINOIJ.N-I) JIMAX0(IB,J-l) DO 3 JMIIB,J XIIX(JM-I) xz-X(JM) X3 - X(JM+I) TERMIIY(JM-I)/((XI-X2)*(XI-X3)) TERMz-Y(JM)/((x2-XI)*(x2-X3)) TERM3-Y(JM+I)/((x3-XI)*(x3-X2)) A-TERMI+TERM2+TERM3 BI-(X2+X3)*TERMI-(XI+X3)*TERM2-(XI+X2)*TERM3 CIX2*X3*TERMI+XI*X3*TERM2+XI*X2*TERM3 IF (JM .GT. IB)GO TO 1A CA-A CB-B cc-C GO TO 15 CAIO.S*(A+CA) CB-0.5*(B+CB) cc-o.5*(c+cc) DIFF - x2-SYL . SUM - SUM + CA*((SYL**2.0)*DIFF + SYL*(DIFF**2.0) + + (DIFF**3.0)/3.0) + + CB*(SYL*DIFF + (DIFF**2.0)/2.0) + CC*DIFF CA-A CB-B cc-C SYL-Xz DIFF - XUP-SYL AVINT - SUM + CA*((SYL**2.0)*DIFF + SYL*(DIFF**2.0) + + (DIFF**3.0)/3.0) + + CB*(SYL*DIFF + (DIFF**2.0)/2.0) + CC*DIFF AVINT-AVINT IF (IND .EQ. 1) RETURN IND-I SYLIXUP XUP-XLO XLo-SYL AVINT--AVINT RETURN END SUBROUTINE IC (KEY,Y,R) I60 I70 168 DIMENSION Y(300) INTEGER R.RR.RRR.RSTOP.RSTART COMMON/ICOND/SO.DO,E0 L - MOD(R.2) RRI3*R RSTOP I (RR/2)-2 RSTART I RSTOP + 3 IF (L .NE. O)RSTART I RSTOP + 2 RRRIRR-Z DO I60 lIl.RSTOP,3 Y(I)ISO Y(I+1)-Do Y(I+2)-E0 CONTINUE DO I70 IIRSTART,RRR,3 Y(I)--S0 Y(I+I)-DO Y(I+2)IEO CONTINUE IF (L .EQ. O)RETURN M.- (RR+1)/2 Y(M-I) I 0.0 Y(M) - 0.0 Y(M+l) I 0.0 RETURN END FUNCTION ZIJ(S,D) C0MM0N/UNITS/SM,OF,GDE.EPM,A,ZP.ZM,CON,TM,DIM,CREF,ZREF COMMON/GREEK/c(2A) REAL IC,M,ICP.MP,ICD.MD ICI(ZP*ZM/2.O)*(((ZP+ZM)*(D+(DF/SM))) + (ZM-ZP)*(S+I.O) + + ((ZREF*ZREF*CREF)/(ZP*ZM*SM))) IF (IC.LE.0.0) ICII.OE-IO Y-O+(DF/SM)-S-1.0 z- D + s + 1.0 + (OF/SM) IF (Y.GE.0.0) YI-l.0E-IO IF (Z.LE.0.0) z-I.0E-10 SUM-((ZM*ZM)*Y - (ZP*ZP)*Z)/(ZM-ZP) ZPPIC(I)*(ZM*(Y)) + C(2)*(ZM*(Y))**I.5 ZMMIC(3)*(ZP*(Z)) + C(A)*(ZP*(Z))**I.5 ZPMIC(5)*(SUM**I.5) + C(6)*(SUM**2.0) HIC(7) - C(8)/((IC)**0.5) B-C(9) + C(IO)*IC CAI((ZP*ZPP) - (ZM*ZPM))/(Y) CBI((ZP*ZPM) - (ZM*ZMM))/(Z) CE-1.0 - ((0.25)*(ZP/ZM)*((ZP-ZM)**2.0)*ZM*(Y)*(M+B)) CG-1.0 - ((0.25)*(ZM/ZP)*((ZP-ZM)**2.0)*ZP*(Z)*(M-B)) CF-I.0 - ((0.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(Y)*(M+B)) 169 CH-I.0 + ((0.25)*(ZM/ZP)*((ZP*ZP)-(ZM*ZM))*ZP*(Z)*(M-B)) CCI((ZP*ZPP) + (ZM*ZPM))/(Y) CDI((ZP*ZPM) + (ZM*ZMM))/(Z) SUMPI(ZP- ZM)/2. 0 ICPIZP*ZM*(ZM- ZP)/2. 0 ZPPPI-(C(I)*ZM) - (I. 5*ZM*C(2)*(ZM*(Y))**O. 5) ZMMP-+(C(3)*ZP) + (I. 5*ZP*C(A)*(ZP*(Z))**O. 5) ZPMPII.5*C(5)*(SUM**O.5)*SUMP + 2.0*C(6)*(SUM)*SUMP MPI(C(8)/(IC**I.5))*ICP BP'C(IO)*ICP CCPI((ZP*ZPPP) + (ZM*ZPMP))/(Y) + +((ZP*ZPP) + (ZM*ZPM))/((Y)**2.0) CDPI((ZP*ZPMP) + (ZM*ZMMP))/(Z) + -((ZP*ZPM) + (ZM*ZMM))/((Z)**2.0) CFPI-((O.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(Y)*(MP+BP)) + + ( (o . 25) * (ZP/ZM) * ( (ZPAZP) - (ZM*ZM) ) *ZM* (PH-8)) CHPI+((O.25)*(ZM/ZP)*((ZP*ZP)-(ZM*ZM))*ZP*(Z)*(MP-BP)) + +((O.25)*(ZM/ZP)*((ZP*ZP)-(ZM*ZM))*ZP*(M-B)) CAPI((ZP*ZPPP) - (ZM*ZPMP))/(Y) + +((ZP*ZPP) - (ZM*ZPM))/((Y)**2.0) CBPI((ZP*ZPMP) - (ZM*ZMMP))/(Z) + -((ZP*ZPM) - (ZM*ZMM))/((Z)**2.0) CEPI-((O.25)*(ZP/ZM)*((ZP-ZM)**2.O)*ZM*(Y)*(MP+BP)) + +((0.25)*(ZP/ZM)*((ZP-ZM)**2.0)*ZM*(M+B)) COP--((0.25)#(ZM/ZP)*((ZP-ZM)**2.0)*ZP*(Z)*(MP-BP)) + -((0.25)*(ZM/ZP)*((ZP-ZM)**2.0)*ZP*(M-B)) DCLI((ZP*ZP*ZPP) + (2.0*ZM*ZP*ZPM) + (ZM*ZM*ZMM)) SCLI((ZP*ZP*ZPP) - (ZM*ZM*ZMM)) SCLPI(ZP*ZP*ZPPP) - (ZM*ZM*ZMMP) DCLPI((ZP*ZP*ZPPP) + (2.0*ZP*ZM*ZPMP) + (ZM*ZM*ZMMP)) SUMDI-(ZP+ZM)*0.5 lCDIZP*ZM*(ZP+ZM)/2.O ZPPDI(C(I)*ZM) + (I.5*ZM*C(2)*(ZM*(Y))**O.5) ZMMDI(C(3)*ZP) + (1. 5*ZP*C(A)*(ZP*(Z))**O. 5) ZPMDII .5*C(5)*(SUM**0. 5)*SUMD + 2. 0*C(.6)*(SUM)*SUMD MDI(O. 5*C(8)/(IC**I. 5))*|CD BDIC(IO)*ICD CADI((ZP*ZPPD)/(Y)) - ((ZM*ZPMD)/(Y)) + '(((ZP*ZPP) - (ZM*ZPM))/((Y)**2.0)) CBDI((ZP*ZPMD)/(Z)) - ((ZM*ZMMD)/(Z)) + -(((ZP*ZPM) - (ZM*ZMM))/((Z)**2.0)) CEDI-(0.25)*(ZP/ZM)*((ZP-ZM)**2.O)*ZM*(M+B) + - ((0.25)*(ZP/ZM)*((ZP-ZM)**2.0)*ZM*(Y)*(MD+BD)) CGDI-(0.25)*(ZM/ZP)*((ZP-ZM)**2.O)*ZP*(M-B) + - ((O.25)*(ZM/ZP)*((ZP-ZM)**2.0)*ZP*(Z)*(MD-BD)) CCDI((ZP*ZPPD)/(Y)) + ((ZM*ZPMD)/(Y)) + -(((ZP*ZPP) + (ZMAZPM))/((Y)**2.0)) CDDI((ZP*ZPMD)/(Z)) + ((ZM*ZMMD)/(Z)) + -(((ZP*ZPM) + (ZM*ZMM))/((Z)**2.0)) CFDI-(O.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(M+B) + - ((O.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(Y)*(MD+BD)) CHDI-(O.25)*(ZM/ZP)*((ZP*ZP)-(ZM*ZM))*ZP*(M-B) + + ((0.25)*(ZM/ZP)*((ZP*ZP)-(ZM*ZM))*ZP*(Z)*(MD-BD)) 222 170 SCLDI(ZP*ZP*ZPPD) - (ZM*ZM*ZMMD) DCLDI((ZP*ZP*ZPPD) + (2.0*ZP*ZM*ZPMD) + (ZM*ZM*ZMMD)) DO 222 I-A.10 TEST - 0.0 IF (C(I) .NE. 0.0)TEST - 1.0 CONTINUE ZIJI((CA*CE) — (CB*CG))/(ZP*ZM) RETURN ENTRY zss ZIJI((CA*CE) - (CB*CG))/(ZP*ZM) RETURN ENTRY ZDD ZIJI((CF*CC) + (CD*CH))/(ZP*ZM) RETURN ENTRY ZSD ZIJI((-CA*CF) - (CB*CH))/(ZM*ZP) RETURN ENTRY ZDS ZIJI((-CC*CE) + (CD*CG))/(ZM*ZP) RETURN ENTRY ZOE ZIJI-(96A85.0*SM*A*A*DCL)/(ZP*ZM*CON*(8.3IAE7)*TM) RETURN ENTRY ZSE ZIJI(96A85.O*SM*A*A*SCL)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY ZSSP \ ZIJI((CAP*CE) + (CEP*CA))/(ZP*ZM) + -((CB*CGP) + (CG*CBP))/(ZP*ZM) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJ = 0.0 RETURN ENTRY ZODP ZIJI((CF*CCP) + (CC*CFP))/(ZP*ZM) + + ((CD*CHP) + (CH*CDP))/(ZP*ZM) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJ - 0.0 RETURN ENTRY ZSDP ZIJI-((CA*CFP) + (CF*CAP))/(ZP*ZM) + - ((CB*CHP) + (CBP*CH))/(ZP*ZM) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJ - 0.0 RETURN ENTRY ZDSP ZIJI((CD*CGP) + (CG*CDP))/(ZP*ZM) + '((CE*CCP) + (CEP*CC))/(ZP*ZM) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJ - 0.0 RETURN ENTRY ZSEP ZIJI(96A85.O*SM*A*A*SCLP)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY ZOEP . ZIJI-(96A85.O*SM*A*A*DCLP)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY ZSSD 171 ZIJI((CA*CED) + (CE*CAD))/(ZM*ZP) + -((CB*CGO) + (CGICBD))/(ZM*ZPI IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJ = 0.0 RETURN ENTRY ZDDD ZIJI((CF*CCD) + (CC*CFD))/(ZM*ZP) + + ((CD*CHD) + (CH*CDD))/(ZM*ZP) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJ - 0.0 RETURN ENTRY ZSDD ZIJI-((CA*CFD) + (CF*CAD))/(ZM*ZP) + -((CB*CHD) + (CBD*CH))/(ZM*ZP) ' IF (c(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJ = 0.0 RETURN ENTRY ZDSD - ZIJI((CD*CGD) + (CDDICG))/(ZM*ZP) + -((CC*CED) + (CE*CCD))/(ZM*ZP) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJ - 0.0 RETURN ENTRY ZSED ZIJI(96A85.O*SM*A*A*SCLD)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY ZOED ZIJI-(96A85.O*SM*A*A*DCLD)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN END FUNCTION ZIJDN(S.D) C0MM0N/UNITS/SM,OF,GDE.EPM,A.ZP.ZM.CON,TM,DIM,CREF,ZREF C0MM0N/GREEK/C(2A) REAL IC,M,ICP,MP,ICO,MD ICI(ZP*ZM/2.O)*(((ZP+ZM)*(D+(DF/SM))) + (ZM-ZP)*(S+I.O) + + ((ZREF*ZREF*CREF)/(ZP*ZM*SM))) IF (IC.LE.0.0) ICIl.0E-IO Y-D+(DF/SM)-S-I.o z- D + S + 1.0 + (OF/SM) IF (Y.GE.0.0) Y--l.OE-IO IF (Z.LE.0.0) z-1.0E-10 SUMI((ZM*ZM)*Y - (ZP*ZP)*Z)/(ZM-ZP) ZPPIC(I)*(ZM*(Y)) + C(2)*(ZM*(Y))**I.5 ZMMIC(3)*(ZP*(Z)) + C(A)*(ZP*(Z))**I.5 ZPMIC(5)*(SUM**I.5) + C(6)*(SUM**2.0) M-C(7) - C(8)/((IC)**0.5) B-C(9) + C(IO)*|C CAI((ZP*ZPP) - (ZM*ZPM))/(Y) CBI((ZP*ZPM) - (ZM*ZMM))/(Z) CE-1.0 - ((0.25)*(ZP/ZM)*((ZP-ZM)**2.0)*ZM*(Y)*(M+B)) CGII.O - ((0.25)*(ZM/ZP)*((ZP-ZM)**2.0)*ZP*(Z)*(M-B)) CFII.O - ((0.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(Y)*(M+B)) 172 CH-I.0 + ((0.25)*(ZM/ZP)*((ZP*ZP)-(ZMRZM))*ZP*(Z)*(M-B)) CCI((ZP*ZPP) + (ZM*ZPM))/(YI CDI((ZP*ZPM) + (ZM*ZMM))/(Z) SUMPI(ZP-ZM)/2.0 ICPIZP*ZM*(ZM-ZP)/2.O ZPPP--(C(1)*ZM) - (I.5*ZM*C(2)*(ZM*(Y))**O.5) ZMMPI+(C(3)*ZP) + (l.5*ZP*C(A)*(ZP*(Z))**0.5) ZPMPII.5*C(5)*(SUM**O.5)*SUMP + 2.0*C(6)*(SUM)*SUMP MP'(C(3)/(IC**I.5))*ICP BPIC(IO)*ICP CCPI((ZP*ZPPP) + (ZMAZPMP))/(Y) + +((ZP*ZPP) + (ZM*ZPM))/((Y)**2.0) CDPI((ZP*ZPMP) + (ZM*ZMMP))/(Z) + -((ZP*ZPM) + (ZM*ZMM))/((Z)**2.0) CFPI-((0.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(Y)*(MP+BP)) + +((0.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(M+B)) CHP-+((0.25)*(ZM/ZP)*((ZP*ZP)-(ZMAZM))AZPA(Z)*(MP-BP)) + +((O.25)*(ZM/ZP)*((ZP*ZP)-(ZM*ZM))*ZP*(M-B)) CAPI((ZP*ZPPP) - (ZM*ZPMP))/(Y) + +((ZP*ZPP) - (ZM*ZPM))/((Y)**2.0) CBPI((ZP*ZPMP) - (ZM*ZMMP))/(Z) ‘ + -((ZP*ZPM) - (ZM*ZMM))/((Z)**2.0) CEPI-((O.25)*(ZP/ZM)*((ZP-ZM)**2.O)*ZM*(Y)*(MP+BP)) + +((O.25)*(ZP/ZM)*((ZP-ZM)**2.0)*ZM*(M+D)) COPI-((0.25)*(ZM/ZP)*((ZP-ZM)**2.O)*ZP*(Z)*(MP-BP)) + -((0. 25)*(ZM/ZP)*((ZP- ZM)**2. 0)*ZP*(M- B)) DCLI((ZP*ZP*ZPP) + (2. OAZMIZPAZPM) + (ZM*ZM*ZMM)) SCLI((ZP*ZP*ZPP) - (ZM*2M*ZMM)) SCLPI(ZP*ZP*ZPPP) - (ZM*ZM*ZMMP) DCLPI((ZP*ZP*ZPPP) + (2.0*ZP*ZM*ZPMP) + (ZM*ZM*ZMMP)) SUMDI-(ZP+ZM)*O.5 ICDIZP*2M*(ZP+ZM)/2.O ZPPDI(C(I)*ZM) + (I.5*ZM*C(2)*(ZM*(Y))**O.5) ZMMDI(C(3)*ZP) + (I.5*ZP*C(A)*(ZP*(Z))**O.5) ZPMDII.5*C(5)*(SUM**O.5)*SUMD + 2.0*C(6)*(SUM)*SUMD MDI(O.5*C(8)/(IC**I.5))*ICD BDIC(IO)*ICD CAD-((ZP*ZPPD)/(Y)) - ((ZM*ZPMD)/(Y)) + '(((ZP*ZPP) - (ZM*ZPM))/((Y)**2.0)). CBD-((ZP*ZPMD)/(Z)) - ((ZM*ZMMD)/(z)) + -(((ZP*ZPM) - (ZH*ZHM))/((Z)**2.0)) CED--(0.25)*(ZP/ZM)*((ZP-ZM)**2.0)*ZM*(M+B) + - ((0.25)*(ZP/ZM)*((ZP-ZM)**2.0)*ZM*(Y)*(MD+BD)) CGDI-(0.25)*(ZM/ZP)*((ZP-ZM)**2.O)*ZP*(M-B) + - ((O.25)*(ZM/ZP)*((ZP-ZM)**2.0)*ZP*(Z)*(MD-BD)) CCDI((ZP*ZPPD)/(Y)) + ((ZM*ZPMD)/(Y)) + -(((ZP*ZPP) + (ZM*ZPM))/((Y)**2.0)) CDDI((ZP*ZPMD)/(Z)) + ((ZM*ZMMD)/(Z)) + -(((ZP*ZPM) + (ZM*ZMH))/((Z)**2.0)) CFDI-(O.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(M+B) + - ((0.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(Y)*(MD+BD)) CHDI-(O.25)*(ZM/ZP)*((ZP*ZP)-(ZM*ZM))*ZP*(M-B) + + ((O.25)*(ZM/ZP)*((ZP*ZP)-(ZM*ZM))*ZP*(Z)*(MD-BD)) 222 173 SCLDI(ZP*ZP*ZPPD) - (ZM*ZM*ZMMD) DCLDI((ZP*ZP*ZPPD) + (2.0*ZP*ZM*ZPMD) + (ZM*ZM*ZMMD)) DO 222 I-A,10 TEST I 0.0 IF (C(I) .NE. 0.0)TEST I I.0 CONTINUE ZIJDNI((CA*CE) - (CB*CG))/(ZP*ZM) RETURN ENTRY ZSSDN ZIJDNI((CA*CE) - (CB*CG))/(ZP*ZM) RETURN ENTRY ZDODN ZIJDNI((CF*CC) + (CD*CH))/(ZP*ZM) RETURN ENTRY ZSDON ZIJDNI((-CA*CF) - (CB*CH))/(ZM*ZP) RETURN ENTRY ZDSDN ZIJDNI((-CC*CE) + (CD*CG))/(ZM*ZP) RETURN ENTRY ZDEDN ZIJDNI-(96A85.O*SM*A*A*DCL)/(ZP*ZM*CON*(8.3IAE7)*TM) RETURN ENTRY ZSEDN ZIJDNI(96A85.0*SM*A*A*SCL)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY ZSSPDN ZIJDNI((CAP*CE) + (CEPICAII/(ZP*ZM) + -((CB*CGP) + (CG*CBP))/(ZP*ZM) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJDN - 0.0 RETURN ENTRY ZDOPDN ZIJDNI((CF*CCP) + (CC*CFP))/(ZP*ZM) + + ((CD*CHP) + (CH*CDP))/(ZP*ZM) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJDN - 0.0 RETURN ENTRY ZSDPON ZIJDNI-((CA*CFP) + (CF*CAP))/(ZP*ZM) + - ((CB*CHP) + (CBP*CH))/(ZP*ZM) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJDN - 0.0 RETURN ENTRY ZDSPON . ZIJDNI((CD*CGP) + (CG*CDP))/(ZP*ZM) + -((CE*CCP) + (CEP*CC))/(ZP*ZM) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJDN - 0.0 RETURN ENTRY ZSEPDN ZIJDNI(96h85.0*SM*A*A*SCLP)/(C0N*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY ZDEPDN ZIJDNI-(96h85.0*SM*A*A*DCLP)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY ZSSDDN 174 ZIJDNI((CA*CED) + (CE*CAD))/(ZM*ZP) + -((CB*CGD) + (CG*CBD))/(ZM*ZP) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJDN - 0.0 RETURN ENTRY ZDDDDN ZIJDNI((CF*CCD) + (CC*CFD))/(ZM*ZP) + + ((CD*CHD) + (CH*CDD))/(ZM*ZP) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJDN - 0.0 RETURN ENTRY ZSDDDN ZIJDNI-((CA*CFD) + (CF*CAD))/(ZM*ZP) + -((CB*CH0) + (CBD*CH))/(ZM*ZP) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJDN = 0.0 RETURN ENTRY ZDSOON ZIJDNI((CD*CGD) + (CDD*CG))/(ZM*ZP) + -((CC*CED) + (CEICCD))/(ZM*ZP) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJDN = 0.0 RETURN ENTRY ZSEDDN ZIJDNI(96A85.O*SM*A*A*SCLD)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY ZDEDDN ZIJDNI-(96A85.0*SM*A*A*DCLD)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN END FUNCTION ZIJUP(S.D) C0MM0N/UNITS/SM,DF.GOE.EPM,A.ZP,ZM.CON,TM,DIM.CREF,ZREF COMMON/GREEK/C(2A) REAL IC,M,ICP,MP,ICD,MO lCI(ZP*ZM/2.O)*(((ZP+ZM)*(D+(DF/SM))) + (ZM-ZP)*(S+I.0) + + ((ZREF*ZREF*CREF)/(ZP*ZM*SM))) IF (IC.LE.0.0) ICII.OE-IO Y-D+(OF/SM)-s-I.0 ZI D + S +1.0 + (DP/SM) IF (Y.GE.0.0) Y--I.0E-IO IF (z.LE.0.0) z-I.0E-10 SUMI((ZM*ZM)*Y - (ZP*ZP)*Z)/(ZM-ZP) ZPPIC(I)*(ZM&(Y)) + C(2)*(ZM*(Y))**I.5 ZMMIC(3)*(ZP*(Z)) + C(A)*(ZP*(Z))**I.5 ZPMIC(5)*(SUM**I.5) + C(6)*(SUM**2.0) MIC(7) - C(8)/((IC)**0.5) B-C(9) + C(IO)*IC CAI((ZP*ZPP) - (ZM*ZPM))/(Y) CDI((ZP*ZPM) - (ZM*ZMM))/(Z) CE-1.0 - ((0.25)*(ZP/ZM)*((ZP-ZM)**2.0)*ZM*(Y)*(M+B)) CG-I.0 - ((0.25)*(ZM/ZP)*((ZP-ZM)**2.0)*ZP*(Z)*(M-B)) CF-1.0 - ((0.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(Y)*(M+B)) 175 CH-I.0 + ((0.25)*(ZM/ZP)*((ZPAZP)-(ZMAZM))*ZPA(Z)*(M-B))~ CCI((ZP*ZPP) + (ZM*ZPM))/(Y) CDI((ZP*ZPM) + (ZM*ZMM))/(Z) SUMP-(ZP-ZM)/2.0 ICPIZP*ZM*(ZM-ZP)/2.O ZPPPI-(C(I)*ZM) - (I.5*ZM*C(2)*(ZM*(Y))**0.5) ZMMP-+(C(3)*ZP) + (I.5*ZP*C(A)*(ZP*(Z))**O.5) ZPMPII.5*C(5)*(SUM**O.5)*SUMP + 2.0*C(6)*(SUM)*SUMP MPI(C(8)/(lC**I.5))*ICP BPIC(I0)*ICP CCPI((ZP*ZPPP) + (ZM*ZPMP))/(Y) + +((ZP*ZPP) + (ZM*ZPM))/((Y)**2.0) CDPI((ZP*ZPMP) + (ZM*ZMMP))/(Z) + -((ZP*ZPM) + (ZM*ZMM))/((Z)**2.0) CFPI-((O.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(Y)*(MP+BP)) + +((0.25)*(ZP/ZM)*((ZP*ZP)'(ZM*ZM))*ZM*(M+B)) CHPI+((O.25)*(ZM/ZP)*((2P*ZP)-(ZMIZM))*ZP*(Z)*(MP-BP)) + +((O.25)*(ZM/ZP)*((ZP*ZP)-(ZM*ZM))*ZP*(M-B)) CAPI((ZP*ZPPP) - (ZM*ZPMP))/(Y) + +((ZP*ZPP) - (ZM*ZPM))/((Y)**2.0) CBPI((ZP*ZPMP) - (ZM*ZMMP))/(Z) + -((ZP*ZPM) - (ZM*ZMM))/((Z)**2.0) CEPI-((0.25)*(ZP/ZM)*((ZP-ZM)**2.O)*ZM*(Y)*(MP+BP)) + +((0.25)*(ZP/ZM)*((ZP-ZM)**2.0)*ZM*(M+B)) CGPI-((0.25)*(ZM/ZP)*((ZP-ZM)**2.O)*ZP*(Z)*(MP-BP)) + -((0.25)*(ZM/ZP)*((ZP-ZM)**2.0)*ZP*(M-B)) DCLI((ZP*ZP*ZPP) + (2.0*ZM*ZP*ZPM) + (ZM*ZM*ZMM)) SCLI((ZP*ZP*ZPP) - (ZMIZM*ZMM)) SCLPI(ZP*ZP*ZPPP) - (ZM*ZM*ZMMP) DCLPI((ZP*ZP*ZPPP) + (2. O*ZP*ZM*ZPMP) + (ZM*ZM*ZMMP)) SUMDI-(ZP+ZM)*O. 5 ICDIZP*ZM*(ZP+ZM)/2.0. ZPPDI(C(I)*ZM) + (l.5*ZM*C(2)*(ZM*(Y))**O.5) ZMMDI(C(3)*ZP) + (l.5*ZP*C(A)*(ZP*(Z))**O.5) ZPMDII.5*C(5)*(SUM**O.5)*SUMD + 2.0*C(6)*(SUM)*SUMD MDI(O.5*C(8)/(IC**I.5))*ICD BDIC(IO)*ICD CAD-((ZP*ZPPD)/(Y)) - ((ZM*ZPMD)/(Y)) + ‘(((ZP*ZPP) - (ZM*ZPM))/((Y)**2.0)) CBDI((ZP*ZPMD)/(Z)) - ((ZMIZMMD)/(Z)) + -(((ZP*zPM) - (ZM*ZMM))/((Z)**2.0)) CED--(O.25)*(ZP/ZM)*((ZP-ZM)**2.0)*ZM*(M+D) + - ((O.25)*(ZP/ZM)*((ZP-ZM)**2.0)*ZM*(Y)*(MD+BD)) CGDI-(O.25)*(ZM/ZP)*((ZP-ZM)**2.O)*ZP*(M-B) + - ((0.25)*(ZM/ZP)*((ZP-ZM)**2.0)*ZP*(Z)*(MD-BD)) CCDI((ZP*ZPPD)/(Y)) + ((ZMAZPMD)/(Y)) + -(((ZP*ZPP) + (ZM*ZPM))/((Y)**2.0)) CDDI((ZP*ZPMD)/(Z)) + ((ZM*ZMMD)/(Z)) + -(((ZP*ZPM) + (ZM*ZMM))/((Z)**2.0)) CFDI-(O.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(M+B) + - ((O.25)*(ZP/ZM)*((ZP*ZP)-(ZM*ZM))*ZM*(Y)*(MD+BD)) CHDI-(O.25)*(ZM/ZP)*((ZP*ZP)-(ZM*ZM))*ZP*(M-B) + + ((0.25)*(ZM/ZP)*((ZP*ZP)'(ZM*ZM))*ZP*(Z)*(MD-BD)) 222 176 SCLDI(ZP*ZP*ZPPD) - (ZM*ZM*ZMMD) DCLDI((ZP*ZP*ZPPD) + (2.0*ZP*ZM*ZPMD) + (ZM*ZM*ZMMD)) DO 222 I-A,IO TEST - 0.0 IF (C(I) .NE. 0.0)TEST - 1.0 CONTINUE ZIJUPI(96A85.O*SM*A*A*SCL)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY ZSSUP ZIJUPI((CA*CE) - (CB*CG))/(ZP*ZM) RETURN ENTRY ZODUP ZIJUPI((CF*CC) + (CD*CH))/(ZP*ZM) RETURN ENTRY ZSDUP ZIJUPI((-CA*CF) - (CB*CH))/(ZM*ZP) RETURN ENTRY ZDSUP ZIJUPI((-CC*CE) + (CD*CG))/(ZM*ZP) RETURN ENTRY zDEUP ZIJUPI-(96485.O*SM*A*A*DCL)/(ZP*ZM*CON*(8.3IAE7)*TM) RETURN ENTRY ZSEUP ZIJUPI(96A85.0*SM*A*A*SCL)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY ZSSPUP ZIJUPI((CAP*CE) + (CEP*CA))/(ZP*ZM) + '((CB*CGP) + (CG*CBP))/(ZP*ZM) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJUP . 0.0 RETURN ENTRY ZDDPUP ZIJUPI((CF*CCP) + (CC*CFP))/(ZP*ZM) + + ((CD*CHP) + (CH*CDP))/IZP*ZM) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJUP = 0.0 RETURN ENTRY ZSDPUP ZIJUPI-((CA*CFP) + (CF*CAP))/(ZP*ZM) + - ((CB*CHP) + (CBP*CH))/(ZP*ZM) 1F (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJUP - 0.0 RETURN ENTRY ZDSPUP ZIJUPI((CD*CGP) + (CG*CDP))/(ZP*ZM) + ’((CE*CCP) + (CEP*CC))/(ZP*ZM) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJUP - 0.0 RETURN ENTRY ZSEPUP ZIJUPI(96h85.0*SM*A*A*SCLP)/(CON*(8.3I4E7)*TM*ZP*ZM) RETURN ENTRY ZDEPUP ZIJUPI-(96h35.0*SM*A*A*DCLP)/(CON*(8.3I4E7)*TM*ZP*ZM) RETURN ENTRY ZSSDUP 177 ZIJUPI((CA*CED) + (CE*CAD))/(ZM*ZP) + '((CB*CGD) + (CG*CBD))/(ZM*ZP) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJUP - 0.0 RETURN ENTRY ZDDDUP ZIJUPI((CF*CCD) + (CC*CFD))/(ZM*ZP) + + ((CD*CHD) + (CH*CDD))/(ZM*ZP) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJUP - 0.0 RETURN ENTRY ZSDDUP ZIJUPI-((CA*CFD) + (CF*CAD))/(ZM*ZP) + -((CB*CHD) + (CBD*CH))/(ZM*ZP) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJUP - 0.0 RETURN ENTRY ZDSDUP ZIJUPI((CD*CGD) + (CDD*CG))/(ZM*ZP) + -((CC*CED) + (CE*CCD))/(ZM*ZP) IF (C(2) .EQ. 0.0 .A. TEST .EQ. 0.0)ZIJUP - 0.0 RETURN ENTRY ZSEDUP ZIJUPI(96A85.0*SM*A*A*SCLD)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY ZDEDUP ZIJUPI-(96A85.O*SM*A*A*DCLD)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN END FUNCTION DDZIJ(S,D) COMMON/UNITS/SM,OF,GDE,EPM,A,ZP,ZM,CON,TM,DIM,CREF,ZREF C0MM0N/GREEK/C(2A) Y-O+(DF/SM)-s-I.0 z- D + S +1.0 + (DF/SM) IF (Y.GE.0.0) Y--I.0E-10 IF (z.LE.O.0) ZII.OE-IO SUMI((ZM*ZM)*Y - (ZP*ZP)*Z)/(ZM-ZP) ZPPIC(I)*(ZM*(Y)) + C(2)*(ZM*(Y))**I.5 INN-C(3)*(ZP*(Z)) + C(A)*(ZP*(Z))**I-S ZPMIC(5)*(SUM**I.5) + C(6)*(SUM**2.0) DCLI((ZP*ZP*ZPP) + (2.0*ZM*ZP*ZPM) + (ZM*ZM*ZMM)) SCLI((ZP*ZP*ZPP) - (ZM*ZM*ZMM)) SUMP-(ZP-ZM)/2.0 ZPPPI-(C(l)*ZM) - (I.5*ZM*C(2)*(ZM*(Y))**O.5) ZMMP-+(C(3)*ZP) + (I.5*ZP*C(A)*(ZP*(Z))**O.S) ZPMPII.5*C(5)*(SUM**O.5)*SUMP + 2.0*C(6)*(SUM)*SUMP DCLPI((ZP*ZP*ZPPP) + (2.0*ZP*ZM*ZPMP) + (ZM*ZM#ZMMP)) SCLPI(ZP*ZP*ZPPP) - (ZM*ZM*ZMMP) SUMDI-(ZP+ZM)*O.5 ZPPDI(C(I)*ZM) + (I.5*ZM*C(2)*(ZM*(Y))**O.5) ZMMDI(C(3)*ZP) + (l.5*ZP*C(A)*(ZP*(Z))**O.5) 178 ZPMDII.5*C(5)*(SUM**O.5)*SUMD + 2.0*C(6)*(SUM)*SUMD DCLDI((ZP*ZP*ZPPD) + (2.0*ZP*ZM*ZPMD) + (ZM*ZM*ZMMD)) SCLDI(ZP*ZP*ZPPD) - (ZM*ZM*ZMMD) DDZIJI(96A85.O*SM*A*A*SCLP)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY DDZDE DDZIJI-(96A85.0*SM*A*A*DCL)/(ZP*ZM*CON*(8.3IAE7)*TM) RETURN ENTRY DOZSE DDZIJI(96A85.0*SM*A*A*SCL)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY DDZSED DDZIJI(96A85.0*SM*A*A*SCLD)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY DDZDED DDZIJI-(96A85.0*SM*A*A*DCLD)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY DDZSEP DDZIJI(96A85.0*SM*A*A*SCLP)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY OOZOEP DDZIJI-(96A85.O*SM*A*A*DCLP)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN END FUNCTION PPZIJ(S.D) , C0MM0N/UNITS/SM,OF,GOE.EPM,A,ZP,ZM,C0N,TM,OIM,CREF,ZREF C0MM0N/GREEK/C(2A) Y-D+(DF/SM)-S-I.o ZI D + S + 1.0 + (DP/SM) IF (Y.GE.0.0) Y--1.0E-10 IF (z.LE.0.0) z-I.0E-IO SUMI((ZM*ZM)*Y - (ZP*ZP)*Z)/(ZM-ZP) ZPP-C(I)*(ZM*(Y)) + C(2)*(ZH*(Y))**1.5 ZMMIC(3)*(ZP*(Z)) + C(A)*(ZP*(Z))**I.5 ZPMIC(5)*(SUM**I.5) + C(6)*(SUM**2.0) DCLI((ZP*ZP*ZPP) + (2.0*ZM*ZP*ZPM) + (ZM*ZM*ZMM)) SCLI((ZP*ZP*ZPP) - (ZM*ZM*ZMM)) SUMP-(ZP-ZM)/2.0 ZPPPI-(C(I)*ZM) - (I.5*ZM*C(2)*(ZM*(Y))**O.5) ZMMP-+(C(3)*ZP) + (I.5*ZP*C(A)*(ZP*(Z))**O.5) ZPMPII.5*C(5)*(SUM**O.5)*SUMP + 2.0*C(6)*(SUM)*SUMP DCLP-((ZP*ZP*ZPPP) + (2.0*ZP*ZM*ZPMP) + (ZM*ZM*ZMMP)) SCLPI(ZP*ZP*ZPPP) - (ZM*ZM*ZMMP) SUMDI-(ZP+ZM)*O.5 ZPPDI(C(I)*ZM) + (I.5*ZM*C(2)*(ZM*(Y))**O.5) ZMMDI(C(3)*ZP) + (I.5*ZP*C(A)*(ZP*(Z))**O.5) ZPMDII.5*C(5)*(SUM**O.5)*SUMD + 2.0*C(6)*(SUM)*SUMD DCLDI((ZP*ZP*ZPPD) + (2.0*ZP*ZM*ZPMD) + (ZM*ZM*ZMMD)) 179 SCLDI(ZP*ZP*ZPPD) - (ZMIZMIZMMD) PPZIJI(96A85.0*SM*A*A*SCL)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY PPZDE ’ PPZIJI-(96A85.O*SM*A*A*DCL)/(ZP*ZM*CON*(8.3IAE7)*TM) RETURN ENTRY PPZSE PPZIJI(96A85.O*SM*A*A*SCL)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY PPZSED . PPZIJI(96A85.0*SM*A*A*SCLD)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY PPZDED PPZIJI-(96A85.O*SM*A*A*DCLD)/(CON*(8.31AE7)*TM*ZP*ZM) RETURN ENTRY PPZSEP PPZIJI(96A85.O*SM*A*A*SCLP)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN ENTRY PPZOEP PPZIJI-(96A85.O*SM*A*A*DCLP)/(CON*(8.3IAE7)*TM*ZP*ZM) RETURN END. WCES REFERENCES Agar, J. N., v n es 1; Elegtroghemistgx and Electroshgmigal Enginegring (Wiley, New York, 1963). Bauerle, J. E., J. Phys. Chem. Solids, 39, 2657 (1969). Bleaney. B. 1., and B. Bleaney. Elegtrigitx 321 Magnetism (Oxford. Clarendon Press, 1957). . Bockris, J. 0'! and S. Srinivasan, Fuel Cells: Their Elegtroghemigggx (NcGraw-Hill, New York, 1969). Borey, B. K., F. H. Horne, Mat. Res. Soc. Symp., 29, 253 (1982). Boyce, I. 8., and R. C. Diprima, Elementary Differential Bguations gag 923535;! yalue Prob1ems (Riley, New York, 1977). Brillouin, L.. J. Phys. Radium.,,;1, 379 (1956). Braunshtein D., D. S. Tannhauser, and I. Riess, J. Electrochem. 800., 128,82 (1981). Brumlev6 T. R., and R. P. Buck, J. Electroanal. Chem.. 29, 1 (1978). Carslav, B. 8., and J. c. Jaeger, Congugtion o_f Heat 33 Solid; (Oxford University Press, London, 1959). Chang, H.-C, and G. Jaffé, J. Chem. Phys., 39, 1071 (1952). Choudhury N. S., and J. I. Patterson, J. Electrochem. Soc.. 118, 1398 (1971). Christy, R. V., J. Chem. Phys., 3 , 1148 (1961). b Chenng, K. Y., B. C. H. Steele, and G. J. Dudley, Fast _2_ Trgngport in Solid; (Elsevier North Holland, Inc.. 1979). 180 181 Cohen, B., and J. l. Cooley, Biophys. J., 5, 145 (1965). Coven, I. A., C. L. Foiles, and D. L. Bdmunds, Fast 'ng Transport in Selig; (Elsevier North Holland, Inc., 1979). Crank, 1., Th; Iathenatig; g; Diffugion (Oxford Univ. Press, London, England, 1956). deGroot, S. 2., L'Effet Soret (Nbrth~Bolland, 1947). deGroot, S. 2., and P. lazur, Noneguilibrigg Thermogzgamigs (North-Holland, Amsterdam, 1969). Dudley, G. J. K. Y. Cheung, and B. C. B. Steele, J. Sol. St. Chem., 32, 269 (1980). Dudley, G. 1., and B. C. H. Steele, I. Sol. St. Chem., ‘2;, 1 (1977). Dudley, 0. 1., and B. c. 11. Steele, J. Nat. 301., g. 1267 (1978). _ Dudley. G. 1., and B. C. H. Steele, J. 801. St. Chem.. 3;, 233 (1980). Dunlap. P. J. and L. J. Gosting, J. Phys. Chem., 53, 86 (1959). Btsell, T. B., and S. N. Flengas, Chem. Rev.,,19, 339 (1970). Euler, 1., and X. Z. Dehmelt, Z. Elektrochem., §;, 1200 (1957). Feldberg, S. V., J. Phys. Chem., 11, 87 (1970). Foley, R. T., J. Electrochem. 800., ;;§, 13C (1969). Fourier, J. B. 1., ThEorie Analytigue g; 3; Chaleur (Paris, 1822). Haase, 2., Thgrmodygamigs 2; Irreversible Progesseg (Addison-Wesley, 182 Booyman, G. 1., B. Boltan, 1r., P. Mazur, and S. R.' deGroot, Physics, ;9, 1095 (1953). Borne, F. E., 1. Chem. Phys., 15, 3069 (1966). Borne, F. E. Physig; g; Superionig Conductors gag Electrode Mgtgria;s, 1. I. Perram, Ed. (Plenum, New York, 1983). Borne, F. E., and T. G. Anderson, 1. Chem. Phys., 53, 2321 (1970). Borne F. E., 1. B. Leckey, and 1. I. Perram, Phygigs g; Superionig Conductors ang Electrgde Materials, 1. I. Perram, Ed. (Plenum, New York, 1983). BOIard, B. E., and A. B. Lidiard, Reports on Progress in Physics, 21, 161 (1964). 1asinski, E.,,gggh,§gg£gy Bagteries (Plenum Press, New York, 1967). Kahn, D., and 1. N. Mayooek, 1. Chem. Phys., 19,4434 (1966). Kiukkola, K. and C. Wagner, J. Electrochem. Soc., ;_1, 379 (1957). Leckey, 1. E., Ph.D. Dissertation, Michigan State University, 1981. Leckey, 1. E., and F. B. Borne, 1. Phys. Chem., §5, 2504 (1981). Ludwig, C., Akad. 'iss. Wien., 29, 539 (1859). Levine, I. N., Physigal Chemistgz (McGrav-Bill, New York, 1978). Longsworth, L. G., 1. Phys. Chem., 6;, 1557 (1957). Marvell, 1. C., Phil. Trans. Roy. Soc. London, ;55, 459 (1865). Miller, D. G., 1. Phys. Chem., ZQ, 2639 (1966). Ohachi, T., and I. Taniguchi, Acts, 22, 747 (1977). 183 Onsager, L., Phys. Rev., 11, 405 (1931). Onsager, L., Phys. Rev., 38, 2265 (1931). Perram, 1. I., 1. Math. Chen.. Q, 37 (1980). Rice, 2. E., Ph.D. Dissertation, Iichigan State University, 1982. Riess, 1., 1. Electrochem. Soc., 128, 2077 (1981). Rosenberg, D. D., lethgdg £2; ‘32; Numerical Solution 2;, Partial Differential figuations (American Elsevier, New York, 1969). Sandifer, 1. 3.. and n. P. Duck, 1. Phys. Chem., 12. 384 (1976). Schmalzried, 3., Solid State Reactions (Academic Press, New York, London, 1974). Sluyters, 1. H. et. al., Rec. Trav. Chim., 12, 1092, 1101 (1960); g;, 100, .525, 535, 553 (1963); §§, 217, 581, 967 (1964); £1, 729, 740, 751. 764 (1965). Smith, G. D. Egggrigsl Solution ‘2; Partial Differential Eguatigns: Finite Difference leghod (Oxford university Press, Oxford, 1978). Snowden,P. N., and 1. C. R. Turner, Trans. Farad. Soc., 56, 1409, 1812 (1960). Soret, C., Arch. Sci. phys. nat. 2, 48 (1879);§, 209 (1880);C. R. Acad. Sci. Paris, 2;, 289 (1880). Sundheim, B. E., and 1. D. Kellner, 1. Phys. Chen., ‘62,. 1204 (1965). Tannhauser, D. S., J. Electrochem. Soc., 125, 1277 (1978). 9. 4191 Tasaka, N., S. Iorita, and I. Nagasava, 1. Phys. Chem., (1965). Tyrell, H. 1. V., Diffusion and Heat Flow in Liguids (Buttervorth and 184 00., London, 1961). 'agner, C., Progr. Sol. St. Chem., 1, l (1972). Iagner, C., Progr. Sol. St. Chem., 10, 3 (1975). Warburg, E., Drude's Annalen der Physik, .61, 493 (1899); liedemann's Annalen der Physik, 6, 125 (1901). Iatanabe, N., and M. A. V. Devanathan, 1. Electrochem. Soc., 11, 615 (1964). 'einberger, H. F., Partial Differential Bguations (Wiley and Sons, New York, 1965). 'eppner, V. and R. A. Huggins, 1. Electrochem. Soc., ;_1, 1569 (1977). 'hittingham, I. 8., 1. 801. St. Chem., 22, 303 (1979).