(186n) Comparison of Various Techniques for Solving Complex Chemical Equilibrium Problems

Computation of the equilibrium composition in complex reactions give rise to the need to solve systems of non-linear algebraic equations (NLE). These systems of equations are notoriously difficult to solve. Common causes of the difficulties include: 1. Presence of discontinuities and/or regions where some of the functions are undefined, 2. While chemical equilibrium problems have only one physically meaningful solution (with only positive and real components) there may be several additional solutions with negative and/or complex components (in case the NLE solver used supports complex domain computation, e.g., Meintjes and Morgan^[1]). Because of the difficulties involved special techniques and programs were developed for solving such problems (for a current review of available packages see for example Leal et al^[2]). However there are often instances where dedicated software packages cannot be used (for example simulation of a complete process requires repetitive chemical equilibrium computations). In such cases general purpose programming languages or numerical software packages need to be used to solve the chemical equilibrium problems, as well. For example Blecic et al^[3] used the programming language Python for this purpose while de Silva et al^[4] and Cutlip and Shacham^[5] used general purpose NLE solvers available in widely used software packages (such as Excel, Polymath and MATLAB). The objective of this comparison study is to help selecting the best problem formulation, initial estimates for the unknowns and solution method in order to ensure convergence to the one physically feasible solution. A brief description of the various options and techniques that were compared follows.

There are two widely used methods for complex chemical equilibrium calculations: the use of equilibrium constants, and the minimization of the Gibbs energy of the reacting system. Using the equilibrium constants approach, the molar amounts at equilibrium are the unknowns and some of the expressions contain products of unknowns raised to various powers (including non-integer). This structure often yields several solutions that may include negative and/or complex molar amounts. In the Gibbs energy minimization approach the molar amounts at equilibrium and Lagrange multipliers are the unknowns. Logarithms of the molar amounts of all the compounds appear in the equations. Consequently, a negative value of a compound, which is obtained during the iterations may stop the computation if the solver used does not support complex number arithmetic. Because of that, software packages that do not support such arithmetic (Python, Excel solver, Polymath) require special approaches to ensure greater than zero molar amounts during the solution. The CONLES (Constrained NLE Solver) program (Shacham^[6]) was used to test iteration schemes that require positive values of some of the unknowns. This program uses a step length restricted Levenberg - Marquardt algorithm to prevent reaching near-zero or negative values of the molar amounts throughout the iterations. MATLAB supports complex arithmetic and its fsolve function was used for solving the NLE's in such an environment. The “trust-region dogleg” algorithm available of the fsolve function was used for solution.

We have found that a critical requirement for reaching the physically feasible solution is the use of initial estimates for some “key” (molar amount) unknowns that do not violate the atom balance equations. These “key” variables are selected based on their expected (substantial) amount. The selection between several molecules containing one common atom is based on their standard Gibbs energy values. Initial estimates for non-key variables can be often obtained directly by substituting the key variable values into some of the equations.

The test problem used for comparison of the various techniques is the “combustion of propane in air”: C₃H₈ + (R/2)*(O₂ + 4N₂) → Products problem, where R = 10 is the relative air to fuel ratio, as presented by Meintjes and Morgan^[1]. The compounds that assumed to be in the equilibrium mixture at 2200 K are: CO₂, H₂O, N₂, CO, H₂, H, OH, O, NO and O₂. The equilibrium constants are specified and the equilibrium composition of the effluent mixture is to be determined.

Using the “equilibrium constant” formulation involves the calculation of square roots of the molar amounts of CO₂, H₂O, N₂, CO and the total number of moles at equilibrium. Meintjes and Morgan, used a special transformation to reduce the number of equations and to eliminate the square root terms, and then used a “continuation” type method to find all the solutions of the system. They identified a total of 10 solutions. One of them is the physically feasible solution (where the molar amounts of all the compounds are positive), three additional real solutions (where some of the molar amounts are negative) and the other are six complex solutions.

We prepared an additional Gibbs energy formulation of that problem, which requires addition of four Lagrange multiplier unknowns and requires calculation of logarithm of the concentrations of the 10 compounds. Several versions of the two problem formulations were prepared: 1. the equations of the “equilibrium constant” formulation were modified so as to replace square roots of the unknowns with integer powers, 2. the equations of the “Gibbs energy” formulation were modified so as to replace the logarithmic terms with exponents, 3. For both formulations the number of implicit equations was reduced to the number of “key variables”, while the additional equations were explicitly solved for the remaining variables.

Upper and lower bounds for the key variables were determined using the atom balance equations and 100 random sets of initial estimates (within the lower and upper bounds) for these variables were generated. Thirteen various formulations of the problem were solved using the fsolve and CONLES functions starting from the 100 randomly generated initial estimates.

The worse performance was observed in the fsolve solution of the “equilibrium constant” formulation, where the square root terms were eliminated and some of the implicit equations were converted to implicit ones. In this case convergence to the physically feasible solution was obtained only from 7 starting points and convergence to a physically infeasible solution (negative concentration of one of the compounds) was obtained from the remaining 93 starting points. The best results (i.e., convergence to the physical solution from the 100 starting points) were obtained by both fsolve and CONLES for the unmodified forms of both the “Equilibrium” and “Gibbs energy” formulations. The complete set of results of the comparison study will be presented in the conference.

The results of this study were validated by solving additional examples. They show that using the preferred problem formulation and programs with robust initial estimates for key variables, increases? substantially the probability of obtaining the physically feasible solutions of complex chemical equilibrium problems.

References

1. Meintjes, K. and Morgan A. P. ACM. Trans. Math. Software, 16(2) 143-151, 1990
2. Leal A. M. M., Blunt M. J. and LaForce, T. C., Geochim. Cosmochim. Acta. 134, 301-322, 2014
3. Blecic J., Harrington, J. and Bowman M. O., ApJ Supplement Series, 225:4, 2016
4. De Silva A. L., Malfatti C. and Muller I. L., Int. J. Hydrogen Energy, 34, 4321-4330, 2009
5. Cutlip, M. B. and Shacham, M. Problem Solving In Chemical and Biochemical Engineering with Polymath, Excel and MATLAB. Prentice-Hall, Upper Saddle River, New-Jersey, 2008
   1. Shacham M., Int. J. Numer. Meth. Eng. 23, 1455-1481,1986