Outline of the Presentation

1. A Systematic Procedure for Analysis and Modification of Chemical Equilibrium Related Equation Sets to Enable Convergence to Physically Meaningful Solutions Mordechai Shacham Ben-Gurion University Beer-Sheva, Israel Neima Brauner Tel-Aviv University Tel-Aviv, Israel

2. Outline of the Presentation • Demonstration of the two commonly used formulations of the complex chemical equilibrium problems and the expected difficulties in solving the resultant NLE (Nonlinear Algebraic Equation) system using general purpose NLE solver and optimization programs • Presentation of the proposed algorithm for determination of the optimal computational sequence and physically feasible initial guesses that enable convergence to the one, physically feasible solution. • Demonstration of the convergence of the modified equation set when starting the iterations from 100 randomly generated initial guesses (within the feasible region)

3. Combustion of Propane (C3H8) in Air (O2 and N2)* • The reaction is C3H8 + (R/2)(O2 + 4N2) → Products where R is the relative air to fuel ratio • The following 10 compounds are assumed in the equilibrium mixture: CO2, H2O, N2, CO, H2 , H, OH, O, NO, O2 • The equilibrium constants at 2200 K are specified • The equilibrium composition of the effluent mixture is to be calculated *Meintjes and Morgan, ACM. Trans. Math. Software, 16(2), 143-151(1990)

4. Combustion of Propane – Atom Balance and Equilibrium Equations (Polymath* Formulation) The sqrt terms may yield error messages or complex solutions if the associated molar amount becomes negative during the iterations *POLYMATH is a product of Polymath Software (www.polymath-software.com )

5. Combustion of Propane – The Only Physically Feasible Solution

6. Combustion of Propane – The 4 Real Solutions Meintjes and Morgan1 used a special “Continuation” type method2 to find all the realand complex solutions. A total of 10 solutions, out of which only one is physically feasible. 1Meintjes and Morgan, ACM. Trans. Math. Software, 16(2), 143-151(1990) 2Morgan A. P., Solving Polynomial Systems Using Continuation for Scientific and Engineering Problems, Prentice-Hall, 1987

7. Catalytic Cracking of Ethane to Form Hydrogen* • Catalytic steam cracking of Ethane to form hydrogen • at temperature T = 1000 K, pressure of P = 1 atm. • Feed - contains 4 moles of H2O per mole of CH4. • The following 9 compounds are assumed in the equilibrium mixture: • CH4, C2H4, C2H2, CO2, CO, O2, H2, H2O, C2H6 • The compound Gibbs Energy values at T =1000K [kcal/gr-mole ]: • CH4 = 4.61, C2H4 = 28.249, C2H2 = 40.604, • CO2 =-94.61, CO = -47.942, O2 = 0, H2 = 0, H2O = -46.03 • C2H6 = 26.13 • The equilibrium composition of the effluent mixture is to be calculated *Balzisher,R. E. Samuels, M. R. and Eliassen, J. D., Chemical Engineering Thermodynamics, Prentice-Hall, Englewood Cliffs, NJ, 1972

8. Atom Balance Equations 3 atoms balance equations should be satisfied:

9. Minimal Gibbs Energy Formulation The solution with this formulation involves the minimization of the total Gibbs energy given by: • ni - number of moles of component i. • c - total number of compounds • R - ideal gas constant • - Gibbs energy of component i at temperature T. • λj,- Lagrange multipliers associated with constrains (gj) due to the 3 atom balance equations. At the minimum point the partial derivatives of F with respect to ni and λj vanish.

10. The Complete Set of Equations (Gibbs Energy Minimization) There are 12 implicit and 1 explicit algebraic equationswith 13unknowns. 10 unknowns are (strictly) constrained to > 0 values. Most solution methods fail ( reaching ≤ 0 valueof some of the unknowns).

11. Difficulties in Solving the Problem in Gibbs Energy Formulation • Initial estimates need to be provided for 12 unknowns. • For some of the unknowns, very small values are expected because of high positive Gibbs energy values (C2H4, C2H2, C2H6), or competing reactions (O2 vs. CO, H2O and CO2). For such compounds it is impossible to come up with sensible initial estimates. • Inadequate initial estimates lead to negative arguments in the logarithmic terms and consequent failure of the solution technique.

12. Changing to Equilibrium Coefficient Formulation 6 reactions of the species involved are considered: The reaction equations are rewritten Definition of Calculation of the equilibrium coeff. Chemical equilibrium formulation αi - stoichiometric coefficient of chemical species Aiin a reaction involving c species. ni - is the mole fraction of component i

13. The Complete Set of Equations (Equilibrium Formulation) • Initial estimates need to be provided for 9 unknowns. • The highly nonlinear terms: (H2/sum)3 and (CH4/sum)2 lead to the existence of several roots: the physically feasible solution - all unknowns have positive and real values, and several physically infeasibleroots - some negative and/or complex unknowns values. • Inadequate initial estimates may lead to negative arguments in the sqrt(O2/sum) term and consequent failure of the solution technique.

14. Attempting to Solve the Systems of Equations • The “safenewt” algorithm of Polymath (www.polymath-software.com) and the “fsolve” function of MATLAB (www.mathworks.com ) were used for solving the various problems. . • The “safenewt” algorithm uses the Newton’s method with line search to ensure reduction of the norm of the function values in every iteration (Press et al1.). The “fsolve” function uses the “trust-region-dogleg” (Powell2) algorithm 1Press, W.H., P.B. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical Recipes, 2nd Ed., Cambridge University Press, Cambridge, 1992. 2Powell, M. J. D., Ch.7 in Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, ed., (1970).

15. Determination of Initial Estimates for the Unknowns 1Balzisher, R. E. Samuels, M. R. and Eliassen, J. D., Chemical Engineering Thermodynamics, Prentice-Hall, Englewood Cliffs, NJ, 1972 2Upper bound is determined on atom balance and Gibbs energy related considerations.

16. Attempting to Solve the Systems of Equations from the Initial Estimate of Balzisher et al. • Polymath (safenewt)–Gibbs Energy formulation. Error message –”attempt to calculate logarithm of a negative number” • Polymath (safenewt) – Equilibrium formulation. Error message–”attempt to calculate square root of a negative number” • MATLAB (fsolve) – Gibbs energy formulation – solution is sought in the complex space. Error message – “No convergence after 1200 function evaluations” • MATLAB (fsolve) – Equilibrium formulation– solution is sought in the complex space. Error message – “No solution found –appears to be locally singular”

17. Algorithm for the Optimal Computational Sequence (1) Identify variables (potential “tear variables”) for which initial estimates based on physical considerations can be obtained. e.g., Species expected to have high concentrations at equilibrium are selected as potential tear variables

18. Algorithm for the Optimal Computational Sequence (2) To determine the optimal computational sequence the Equations • - Equations -Variables’ Boolean Matrix is constructed: rows represent the equations, columns the variables. • Presence of a particular variable in a selected equation: indicated by the entry 1 (one) in the corresponding cell. • The row (equation) with the minimal number of variables is selected to enter the computational sequence, provided that all (except one) variables are potential tear variables. • The selected equation is rewritten in an explicit form: tearvariables appear on the r.h.s., calculated variable on the l.h.s. • The selected equation is entered into the computational sequence: i.e., the associated equation & variables are removed from the boolean matrix.

19. Algorithm for the Optimal Computational Sequence (3) • Explicit Equations- those where only one unknown variable appears, can be solved for that unknown variable. Preferring a formulation where potential singularity causing expressions, (e.g., logarithm, square root ) were removed. • Remaining implicit equations - for the unknown “tear variables” (for which initial estimates has to be assigned). The number of implicit equations must equal to the number of the tear variables.

20. Equation Sequencing in the Equations–Variables Boolean Matrix Tear variables Calculated variable Equation Sequencing, Step 1 Rewrite the equation: f1 = (H2O*CO)*1.3786 - (CO2 *H2) = 0 in the form CO2 = (H2O * CO) * 1.3786 / H2 Potential tear variables are shown in bold letters.

21. Equation Sequencing, Step 2 Tear variable Calculated variable Rewrite the equation f5: f5 = (CO / sum * (H2 / sum) ^ 3) - 26.6291 * CH4 * H2O / sum ^ 2 = 0in the form CH4 = (CO / sum * (H2 / sum) ^ 3) * sum ^ 2/(26.6291 * H2O)

22. Equation Sequencing – Final steps • f4 can be solved for O2 and f6 for C2H2, then f3 can be solved for C2H4 and f2 for C2H6. • To eliminate the possibility of calculating square root of a negative number: f4 = H2O - 1.1473e10 * H2 *sqrt(O2 / sum) Changed to : f4 = H2O ^ 2 - 1.1473e10^2 * H2^2 * (O2 / sum) = 0

23. The Complete Modified Set of Equations (Equilibrium Formulation) • Initial estimates need to be provided only for 4 unknowns. Those estimates can be obtained based on physical considerations. • The sqrt(O2/sum) term has been eliminated: Search of solutions in the complex space is prevented.

24. Equilibrium Composition and Function Values (Equilibrium Formulation, Polymath Solution) The concentration of some of the species is practically zero. During the solution process <=0 values can be attained.

25. The Complete Modified Set of Equations (Gibbs Energy Formulation) • 4 implicit and 9 explicit algebraic equations • Logarithms need to be calculated only for variables for which good initial estimates are available

26. Equilibrium Composition and Function Values (Gibbs Energy Formulation, Polymath Solution) The solution matches that of the “Equilibrium Formulation” (with up to six decimal digits).

27. The Effect of the Initial Estimates on Convergence • To check the sensitivity of the convergence to the initial guess used, 100 random starting points were generated between the minimal and maximal values of the tear variables. • The fsolve function of MATLAB was used to solve the system of equations for both the Gibbs Energy and the Equilibrium formulations. • Using the proposed algorithm, convergence to the correct solution was achieved from the 100 initial guesses within maximum of 10 iterations.

28. Conclusions • Complex chemical equilibrium problems may be practically impossible to solve by nonlinear solvers available in general purpose numerical computation packages (e.g. Polymath and MATLAB), when using the basic theory-based formulation of the equations and arbitrary initial estimates for the unknowns. • The use of a new procedure has been suggested that involves selection of tear variables, equation sequencing and modifications to remove singularities and reduce the number of variables that need to be iterated upon. • It was shown that with the proposed procedure, very fast convergence to the accurate physical solutions can be reached from a wide range of initial estimates for the tear variables.