Journal of ThermodynamicsVolume 2012 (2012), Article ID 870631, 18 pagesdoi:10.1155/2012/870631

Research Article

IGE Model: An Extension of the Ideal Gas Model to Include Chemical Composition as Part of the Equilibrium State

Christopher P. Paolini and Subrata Bhattacharjee

Department of Mechanical Engineering, San Diego State University, San Diego, CA 92182, USA

Received 14 August 2011; Accepted 1 November 2011

Academic Editor: Hugo Segura

Copyright © 2012 Christopher P. Paolini and Subrata Bhattacharjee. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The ideal gas (IG) model is probably the most well-known gas models in engineering thermodynamics. In this paper, we extend the IG model into an ideal gas equilibrium (IGE model) mixture model by incorporating chemical equilibrium calculations as part of the state evaluation. Through a simple graphical interface, users can set the atomic composition of a gas mixture. We have integrated this model into a thermodynamic web portal TEST (http://thermofluids.sdsu.edu/) that contains Java applets for various models for properties of pure substances. In the state panel of the IGE model, the known thermodynamic properties are entered. For a given pressure and temperature, the mixture's Gibbs function is minimized subject to atomic constraints and the equilibrium composition along with thermodynamic properties of the mixture are calculated and displayed. What is unique about this approach is that equilibrium computations are performed in the background, without requiring any major change in the familiar user interface used in other state daemons. Properties calculated by this equilibrium state daemon are compared with results from other established applications such as NASA CEA and STANJAN. Also, two different algorithms, an iterative approach and a direct approach based on minimizing different thermodynamic functions in different situation, are compared.

1. Introduction

The accuracy of a solution to an engineering thermodynamics problem lies in large part to the material model chosen to model a working fluid. Generally speaking, there are five approaches one can take when modeling a working substance: one can choose to model a substance as follows:(i)a material that undergoes a change in phase,(ii)a condensate (i.e., a solid or liquid substance),(iii)an ideal gas,(iv)a perfect gas, (v)a real gas.

The particular model chosen will dictate which set of equations to use when calculating thermodynamic states. Generally, the simpler the model chosen, the simpler the calculation becomes. The simplest of all models is the perfect gas model, and this model is usually selected when performing state approximations with pencil and paper. But, simplicity tends to come at the expanse of accuracy. The more complex models require either tables of precomputed results or a computer program that implements a numerical algorithm to obtain a solution.

Computing the state of substances that undergo a change in phase requires saturation and superheated tables. For example, when modeling a thermodynamic cycle that uses water as a working fluid, Steam Tables are used to determine the saturation pressure, specific volume, internal energy, enthalpy, and entropy of steam at a given temperature. Other substances such as R-12, NH₃, and R-134 a, should be modeled using a phase change model since there is likely the possibility of a phase transformation during a cycle.

A condensate is a working substance that resides in a solid or liquid phase for the duration of a cycle. Condensates can be modeled as materials that have a constant density and a constant heat capacity which is equal at constant pressure or constant volume.

The state of an ideal gas is modeled using the ideal gas equation, In this model, heat capacity is temperature dependent. The perfect gas model obeys the ideal gas equation; however, heat capacity is considered constant, which renders analytical integration over a temperature range easier to perform when approximating a change in enthalpy or entropy. In this sense, a perfect gas can be considered as a simplified ideal gas.

The real gas model uses a generalized compressibility chart where (1) is modified such that with representing a compressibility factor. The real gas model uses the compressibility factor to account for how a gas deviates from the ideal gas model at low temperatures or high pressures. A table of pre-computed compressibility factors is used to determine the proper value of when computing the state of a particular gas.

The ideal gas law was first introduced by Clapeyron in his 1834 work entitled Mémoire sur la puissance motrice de la chaleur “Memoir on the motive power of heat.” This model assumes that a gas consists of infinitesimally tiny particles where the length scale of a gas molecule diameter is much smaller than the length scale of distances traveled when intermolecular collisions occur. The ideal gas model also assumes that gas molecules interact only through brief, infrequent, and elastic collisions [1]. Thus, under this model, when two or more gas molecules collide, the total translational kinetic energy of the molecules is conserved. At low temperatures and high pressures where molecules have a higher probability of interaction, the ideal gas law is less accurate and fails to adequately describe the state of a gaseous species. Many attempts have been made over the years to improve the ideal gas law to account for molecular interaction. Most notable was J. D. van der Waals who, in 1873, modified the ideal gas law by accounting for the effective volume gas molecules have to occupy by subtracting from the volume where is the volume of one mole of molecules. Second, van der Waals added another term to account for molecular attractive forces which varies proportionally to the inverse square of the distance between molecules. With representing the constant of proportionality, van der Waals’ model is given by For a given species, the constants and are determined by making use of the fact that the critical point lying on a critical isotherm on a pressure versus molar volume plot is an inflection point. By setting one will obtain two equations in the two unknowns, and , which can then be computed from knowledge of the species’ critical temperature and pressure, ,  , respectively. One deficiency with the van der Waals model is a lack of accounting for intermolecular repulsive forces. Since 1873, many other equations of state (EOS) have been proposed to account for both attractive and repulsive forces. One of the most recognized is the model proposed in 1948 by Redlich and Kwong [2], The Redlich and Kwong model preserves the first term of the van der Waals model, but makes an adjustment to the second by introducing a dependence on temperature and the effective volume to account for repulsive forces. Species specific values for the two constants and are then found using the same technique used in the van der Waals model. Many other equations of state have been proposed (Peng and Robinson [3], Kerrick and Jacobs [4]) with each model having its own strengths and weaknesses. For example, while the Redlich and Kwong model does a good job representing the state of many gasses at low temperatures and high pressures, the Kerrick and Jacobs EOS does a good job at high temperatures and high pressures but not low temperatures and high pressures. Nevertheless, the standard ideal gas law of Clapeyron is the simplest model and is widely used in most combustion research and textbooks for educational purposes. Although we are extending the IG model, the same methodology can be applied to other equations of state where the equilibrium distribution is computed simultaneously with the thermodynamic state.

The computational method employed in this work is based on the direct minimization of Gibbs free energy at either constant temperature and pressure or constant enthalpy and pressure. The direct free energy minimization algorithm was originally proposed by White et al. in 1958 [5] and further developed by Zeleznik and Gordon throughout the 1960s [6–10]. The algorithm was originally developed for computing the chemical equilibria of ideal systems that characterize rocket-propellant and other combustion scenarios. Assuming an ideal-gas behavior of such systems is not unreasonable, due to the high temperature and atmospheric level pressures involved. The IGE model discussed in this work assumes a single gas phase and one or more pure condensed (solid or liquid) phases. Inclusion of a pure condensed phase is performed after the equilibrium composition is computed for the gas phase. Then, a vapor pressure test is performed to determine if the inclusion of a pure condensed phase lowers the system Gibbs energy [8, 11]. Upon adding a condensate from the list of atomically feasible candidate species, the equilibrium composition is computed again. This iterative process is repeated until all atomically feasible condensates are considered for inclusion in the final equilibrium products mixture.

Because our IGE model and software implementation is intended to be used to solve primarily high-temperature/low-pressure gas phase problems that arise in combustion scenarios, the computation of chemical equilibria for non-ideal systems is not supported. Much work, however, has been done in the area of nonideal equilibrium computation. Such techniques usually employ stoichiometric and extent of reaction equations to solve multiphase equilibria in nonideal systems. Of particular importance is the work of Castier et al. [12] and Michelsen [13, 14] in which a numerical technique is presented that simultaneously computes the equilibrium distribution while discovering all thermodynamically feasible phases, without requiring the user to prespecify the number and type of phases existing in the equilibrium state. Castier’s method employs a stoichiometric formulation and demonstrates a second-order convergence rate. Nonidealities are handled by reformulating reaction extent equations to include a fugacity and activity coefficient model for gas and liquid phase calculations, respectively. Extent of reaction methods used to determine chemical and phase equilibria in nonideal systems has also been explored by Economou et al. [15]. Economou’s method is based on computing the equilibrium constant for a reaction as a function of the fugacity of each component and then using an EOS that accounts for nonideal behavior to compute fugacity coefficients. Economou shows how a component’s fugacity coefficient can be calculated by expressing an EOS in terms of a compression factor that includes an ideal constituent, a contribution that accounts for molecular repulsive forces, and a contribution that accounts for molecular attractive forces. Economou et al. applied their model using the Peng-Robinson EOS as well as the Soave modification to the Redlich-Kwong EOS [16]. The Redlich-Kwong-Soave EOS given by includes an additional factor that is multiplied by the attractive parameter in (5) and is a function of temperature and an acentric factor that characterizes the nonsphericity of molecules. Application of the model shows how both equations of state accurately capture the effect of high pressure and low temperature on reaction extent and phase-equilibria computations. A modification of the Peng-Robinson EOS can be made that allows the model to better capture vapor pressures of pure fluids. The so-called PRSV2 EOS [17] reformulates the attractive term in the Peng-Robinson EOS and is given by where three parameters, , , and are used to compute a factor that modifies the attractive term factor . In this case, is a function of and a reduced temperature with being a function of the acentric factor . Optimal values for , , and for many fluids are found in the literature [18–20]. Llano-Restrepo and Muñoz-Muñoz [21] employed the PRSV2 EOS with the UNIQUAC activity coefficient model [22] and the Wong-Sandler (WS) mixing rules [23] to compute the vapor-liquid equilibrium (VLE) of the ethylene-ethanol-water ternary system. Llano-Restrepo et al. uses the flash method described in the work by Sandler [23] to compute the equilibrium distribution of the liquid and vapor phases of ethylene-ethanol-water mixtures at 200°C and pressures ranging from 30 to 154 atm. Results showed their numerical model produces mole fraction predictions that compare well with experimental data obtained from the hydration of ethylene for synthetic ethanol production. The Beattie-Bridgeman EOS and the Benedict-Webb-Rubin EOS are frequently used in mechanical engineering for modeling gas phase species. The Beattie-Bridgeman EOS [24] is given by where and are empirical parameters, it is fairly accurate for species with densities up to [25]. The Benedict-Webb-Rubin EOS [26], where the constants , , , , , , , and have been tabulated for many gas species [27], it is accurate for species with densities up to . Nonetheless, the ideal gas model is the most straightforward EOS and is routinely employed to solve most gas phase problems encountered in combustion products equilibria. For this reason, this work is based on extending the ideal gas EOS in our online thermodynamics software suite named TEST to compute the thermodynamic state of the products of combustion reactions based on a chemical equilibrium distribution computed assuming an ideal gas EOS.

2. The IG and IGE Models in Test

The Expert System for Thermodynamics, abbreviated TEST, is an Internet portal for accessing thermodynamics web applications [28, 29]. TEST is freely accessible via URL http://www.thermofluids.net/ and combines applications with multimedia educational content, such as example problems with solutions that are typical of those encountered by university engineering students. In addition, one can access a variety of thermodynamic charts and tables and fifteen chapters of animations that illustrate thermodynamic systems and fundamental concepts. TEST Web applications are called daemons and are mostly implemented as Java applets. Recently, a new kind of web application has been introduced into TEST, termed a Rich Internet Application or RIA, which is implemented using Adobe Flash [30]. Because these Web applications are accessible online, they require nothing more than a standard Web browser to operate, which makes them convenient for researchers, educators, and students who can use the daemons to solve fairly complex engineering problems, quickly and efficiently, without having to download and install platform-specific software. TEST’s daemons are specially designed to allow a user to easily specify a thermodynamic problem and then quickly solve follow-up or closely related problems known as what if? scenarios. For example, using the TEST Refrigeration Cycle daemon, one could configure a problem to solve for the cooling power in tons of a refrigeration system that uses the CFC Freon-12 (R-12 or dichlorodifluoromethane) as a working fluid. After specifying known inputs to an initial problem, such as the mass flow rate of the refrigerant and the temperature of the surroundings, and a solution is obtained, one can quickly modify an input variable and obtain a solution to a related problem. For example, the user may ask “What if R-12 was replaced with a more environmentally friendly refrigerant such as R-134 a? How would the cooling power change?” In this case, the user can simply pull down a menu list of refrigerants, select R-134 a, and click a single button to compute a new solution. TEST features numerous daemons that are specifically designed to solve problems in such areas as IC engines, gas and vapor power cycles, refrigeration, HVAC, combustion, chemical equilibrium, and gas dynamics. Furthermore, TEST features general daemons that can be used to evaluate state properties of various working substances [31] and perform energy, entropy, and exergy analysis of generic open and closed thermodynamic systems.

The latest advancements in TEST have been in the area of chemical equilibrium analysis. The TEST chemical thermodynamic module [32–34] is comprised of a set of combustion and equilibrium analysis daemons for solving both open steady systems problems and unsteady process problems. The combustion daemons can be used for the analysis of premixed or nonpremixed reactants modeled using a perfect gas (PG), ideal gas (IG), or a solid/liquid (SL) model. In the TEST framework, the PG model stipulates that a working fluid obeys the ideal gas equation with heat capacity held constant. In the IG model the working fluid still obeys the ideal gas equation but the heat capacity of the working fluid is a function of temperature. A working fluid modeled using the SL model in TEST has a constant density and constant heat capacity with heat capacity at constant pressure equal to heat capacity at constant volume (i.e., ). TEST provides separate daemons for thermodynamic state evaluation using the IG model when the working fluid is a single species, a binary mixture of two species, or a general mixture of more than two species. To illustrate the IG model and contrast it with the proposed IGE model, the trivial problem of computing the specific enthalpy of oxygen gas over a large temperature range is used as a benchmark. The IG state daemon is accessible by navigating to Daemons States System IG-Model. Figure 1 shows a screen capture of the daemon. The majority of TEST’s daemons follow the look and feel of the interface shown in Figure 1. Each daemon is a graphical Web-based specialized calculator designed for solving a specific class of thermodynamic problem. The ideal gas daemon provides several text entry fields in which a user enters what is known about a particular problem, and a calculate button for computing what is unknown. For example, when evaluating the state of a pure gas, specifying any two thermodynamic properties is sufficient to determine all other properties. Thus a user could enter values for the specific volume and temperature of a selected gas, at which time the daemon would use the ideal gas law to calculate all other properties such as pressure, specific internal energy, specific enthalpy, and specific entropy. The green text entry boxes show input values entered by the user while the cyan boxes show the computed output values. In Figure 1 the user has entered and . The notation and refer to the pressure and temperature, respectively, of the first thermodynamic state specified in the daemon. The @State pull-down menu in the upper-left corner of the daemon allows a user to configure multiple states in which different inputs are specified and then compute all the properties of all the states at once using the blue Super Calculate button. The blue line in Figure 2 shows a plot of the specific enthalpy of oxygen gas at a pressure of 1 atmosphere over the temperature range (298.15 K, 6000 K). The enthalpy values are obtained by configuring the daemon at thirteen different states with different temperatures, T1 through T13, and then using the Super Calculate button to compute thirteen different specific enthalpy values, as indicated by the circle symbols in the plot. The specific enthalpy of oxygen gas for temperatures approximately less than 3000 K Kelvin is accurately shown in the plot. In reality, however, the diatomic oxygen molecule will begin to dissociate at higher temperatures and form monatomic O atoms. The degree of dissociation is temperature and pressure dependent. At temperatures greater than 3000 K, a mixture of O₂ molecules, and O atoms will form and this mixture will have different thermodynamic properties than a pure amount of O₂ molecules. As one can see in Figure 2, the specific enthalpy of the true dissociated mixture of oxygen gas, represented by the red line, rises sharply after 3000 K and is remarkably greater than the specific enthalpy of the pure substance at higher temperatures. The specific enthalpy values shown on the red line were computed using the TEST ideal gas equilibrium of IGE model. The IGE model is the same as the IG model with the exception that the equilibrium composition of the working fluid is computed first and then the thermodynamic state is evaluated based on the computed mixture. Any application in which the IG model can be used, the IGE model can also be used to obtain a more accurate thermodynamic representation of the system.

Figure 1: The single species variant of the TEST IG model system state daemon: computing the specific enthalpy of oxygen gas at 1 atm and 3000 K.

Figure 2: The specific enthalpy of oxygen gas at one atmosphere over the temperature range (298.15 K, 6000 K). The IGE model computes a more accurate thermodynamic state since the true mixture composition is determined before state computation.

The degree of dissociation is shown in Figure 3 where the mole fractions of O₂ molecules and O atoms are plotted with respect to temperature. From this figure one can see that at 6000 K a sample of O₂ molecules will completely dissociate into O atoms. We can use the data computed by the IGE model to examine the endothermicity of the dissociation process. As the temperature of a mixture of molecules and O atoms increases, absorption of heat occurs. The absorbed energy will be distributed into translational (kinetic) energy, internal vibrational and rotational energy, as well as electronic and nuclear energy. If enough vibrational energy is present in an O2 molecule to overcome the double covalent bond holding its two oxygen atoms together, atomization will occur and the O₂ molecule will dissociate into two O atoms. The bond dissociation energy D°  for an OO bond which is broken through the reaction is defined as the standard state enthalpy change for the reaction at a specified temperature, The fraction of O atoms present is dependent on temperature and pressure. At a standard pressure of 1 atmosphere, Figure 3 shows a plot of the mole fraction of O₂ and O for a 1 mole initial mixture of O₂ in thermodynamic equilibrium at a given temperature. The equilibrium reaction for this process is given in (11) which indicates at equilibrium, for a particular temperature, the forward and backward reactions occur at the same rate. In other words, the rate at which O₂ bonds are broken to form oxygen atoms is matched by the rate at which two oxygen atoms recombine to form an oxygen molecule. The amount of O₂ molecules and O atoms present in equilibrium at a particular temperature is defined by the partial pressure equilibrium constant for reaction (11). If we let represent the fraction of each mole of O₂ molecules that dissociate into moles of O atoms, leaving 1- moles of O₂, we can express for reaction (11) in terms of the partial pressures of the reactant O₂ and the product O, Figure 4 shows a plot of the log of versus the reciprocal of temperature obtained from the IGE daemon. One can see a linear relationship with a negative slope is present at high temperatures greater than about 1500 K. For temperatures less than 1500 K, log is constant and independent of temperature. The negative slope indicates reaction (11) is endothermic. As the temperature increases, increases. By Le Chatelier’s Principle the equilibrium will shift to the right in (11) to counteract the stress of additional heat. A right shift will produce an increase in O atoms and a corresponding decrease in molecules. As the number of O atoms increases, the partial pressure of O atoms increases and this increases the equilibrium constant .

Figure 3: The mole fractions of O₂ molecules and O atoms in a one mole initial sample of oxygen gas at standard pressure. As temperature increases, O₂ molecules dissociate into O atoms. At 6000 K, a one mole initial sample of oxygen gas will contain 2 moles of O atoms.

Figure 4: Log of the partial pressure equilibrium constant versus the reciprocal of temperature for reaction (11).

Using the van’t Hoff equation, and the slope of the log v. 1/ plot, we can solve for , the standard state molar enthalpy of reaction, and investigate how varies with temperature. Figure 5 shows a plot of versus temperature. The minima of this function corresponds to the point of intersection in Figure 3 where an equal portion of O₂ molecules and O atoms are present.

Figure 5: Standard state molar enthalpy of reaction versus temperature for reaction (11). which implies the reaction is endothermic.

3. Numerical Method

The theory underlying the computation of the distribution of a mixture of ideal gases in equilibrium at a constant temperature and pressure is based on the minimization of the Gibbs free energy function. The interested reader can find a detailed description of the numerical method used by the TEST IGE model to minimize a system’s Gibbs function in Paolini and Bhattacharjee [35], Gordon and McBride [11, 36, 37], and the classic paper by Gordon et al. [38].

To summarize, the Gibbs function is defined as A system will be in equilibrium when the change in Gibbs free energy of the system vanishes to zero.

Differentiating (14) we obtain At a constant temperature and pressure , (15) reduces to The combined first and second laws of thermodynamics requires which tells us the equilibrium state is also one where the Gibbs free energy of a system has reached the smallest possible or minimum value. From (14) and using the definition of the chemical potential for an ideal gas species , Substituting (18) into (14) and dividing by to obtain a dimensionless form, we obtain but where the last term, , reduces to a constant since the system total pressure is given as an input to the problem. Using (20), (19) now becomes The minimum stationary point of (14) will be the vector of species molar values where vanishes. Differentiating in (14) we obtain but from the isothermal, isobaric Gibbs-Duhem equation we know that and so we seek the unique vector such that Solving (24) amounts to solving a nonlinear constrained minimization problem. The method of Lagrange multipliers is used to solve (24) according to an atomic mass constraint The total number of each atom present in the reactant mixture must be the same in the product mixture at equilibrium since mass is conserved. Because there cannot exist a negative number of moles of a species, we are also bound by the positive moles constraint We can express the atomic constraint equations given by (25) and (26) as equality constraints There will be equations of the form (27) for each atom present in the system. We define the Lagrangian function as The solution of (24) is thus the minimum stationary point where the gradient of (28) vanishes, that is, The TEST IGE model uses a Newton-Raphson technique to solve (29). By expanding (29) we can write In order to satisfy (30), all the terms in the square brackets must vanish to 0. We therefore seek values of Lagrange multipliers such that for all . By observing solving (30) amounts to solving a system of equations for each species where the th species equation is given by and additionally, population constraint equations for each element where the th constraint equation is given by and finally one constraint equation governing the total number of moles present in the system, From (33), (34), and (35), we must solve a system of equations in unknowns. The IGE model uses the Newton-Raphson method to iteratively solve this system of equations. To construct our system matrix we let represent the th species equation of (33), or Expanding using (21) we have By letting we will obtain and our Newton-Raphson equation becomes We will have equations of type (40) for each species . For the population constraint equations given by (34), we let represent the th population constraint equation such that Then and our Newton equation for (34) becomes We will have of these equations for each unique atom present in the system. Finally, we accommodate the total system moles equation in a similar fashion by letting represent (35), where and thus the Newton equation for the total system moles constraint becomes We can represent our system of Newton equations in matrix form by combining (40), (43), and (46) to obtain

Many common problems in combustion, gas turbine analysis, and gas dynamics require the chemical composition of a set of reacting species to be determined not at constant temperature, but at constant enthalpy or entropy. For example, the combustion process of a Diesel cycle can be modeled as a constant pressure adiabatic process where the piston moves to maintain a constant pressure. In combustor design one often is interested in determining the adiabatic flame temperature at constant pressure which amounts to the temperature the products of combustion will attain if no energy is lost to the surroundings. In the simplified analysis of flow through a convergent-divergent nozzle, the gas flow is assumed to be isentropic. The IGE model can be employed to solve these types of problems by constraining enthalpy or entropy, respectively, instead of temperature. To solve an adiabatic problem, the IGE model enforces a constant enthalpy constraint such that a distribution is found that minimizes a Gibbs function such that the enthalpy of the product distribution is equal to the enthalpy of the initial reacting mixture.

The TEST IGE model can solve a constant pressure, constant enthalpy (“”) or constant pressure, constant entropy (“”) equilibrium problem in two different ways. The first is a direct method in which an additional constraint equation is incorporated into the iteration matrix given in (47). The second is an iterative approach that repeatedly solves a constant pressure, constant temperature (“”) problem until the specific enthalpy of the product mixture equals the specific enthalpy of the reactants. The direct constraint equation is derived as follows. Let represent the dynamically changing specific enthalpy of the product mixture and represent the specific enthalpy of the reactants, a constant. Upon convergence of (47) we desire The specific enthalpy of the product mixture is determined after each iteration of (47) through the summation where is the standard state specific molar enthalpy of species . As the number of moles of each species varies, the current value of the product mixture enthalpy varies as For an problem, temperature varies and approaches the adiabatic flame temperature as (47) converges. As varies, the product mixture enthalpy also varies as where is the standard state specific molar heat capacity of species . Equation (51) is clear sincefor an ideal gas. The standard state specific molar enthalpy and heat capacity in (50) and (51) are cast into dimensionless form by dividing by and , respectively. Our additional constraint equation is thus or Equation (53) is inserted into (47) and repeatedly solved yielding correction terms and. The final temperature will be the adiabatic flame temperature and the final product mixture will have a specific enthalpy equal to that of the reactants. The procedure is similar for problems. In this case we require the entropy of the product mixture to equal the entropy of the reactants, or where represents the dynamically changing specific entropy of the product mixture and represents the specific entropy of the reactants, a constant. The specific entropy of the product mixture is determined after each iteration of (47) through the summation where is the specific molar enthalpy of species given by and is the standard state specific molar entropy of species given by As the number of moles of each species varies, the current value of the product mixture entropy varies as For an problem, temperature varies and approaches an isentropic temperature as (47) converges. As varies, the product mixture entropy also varies as The last equality in (59) is from (46). Our additional constraint equation for an isentropic, isobaric process is thus The second method used by the TEST IGE model to solve a or problem is an iterative method in which repeated constant pressure, constant temperature (i.e., ) computations are made until a converged value of temperature is obtained. Pseudocode to compute an adiabatic flame temperature using this iterative method is as in Algorithm 1.

Set

Specify the specific enthalpy and pressure of the reactants

Compute the temperature as a function of and

Compute the equilibrium distribution of the product mixture

    according to and

Repeat

    Set

    Compute the thermodynamic state of this product mixture

    based on and giving a new temperature

    Set

    Compute the equilibrium distribution of the product

    mixture according to and

Until

Algorithm 1

where is a convergence error tolerance. NASA CEA [11, 36] uses the direct approach while STANJAN [39] performs repeated iterations.

4. Model Comparison

To illustrate how the TEST IGE model can extend the capabilities of the IG model we examine how the two models can be used to compute adiabatic flame temperature. For combustion applications using the IG model, TEST includes different suites of daemons for open and closed processes. In each suite, IG model daemons exist for problems in which the reacting fuel and oxidizer components are premixed, or problems in which a fuel and oxidizer enter a combustion chamber separately. In both instances, the IG model assumes the product composition is frozen. That is, one must completely define the distribution of products of combustion prior to computation and the adiabatic flame temperature computed is based on this prespecified distribution. Figure 6 shows a screenshot of the TEST Specific, Closed Process, Premixed, Reaction daemon using the IG Mixture Model. This combustion daemon is designed to compute the state of premixed fuel and oxidizer reactants using the ideal gas law. To use this daemon, one first defines a combustion reaction by specifying reactant and product species as well as the reaction stoichiometry. One can balance, scale, and normalize a reaction using the commands shown in a drop-down menu. After defining a combustion reaction, one navigates to the State Panel where a sufficient number of known thermodynamic properties needed to compute the thermodynamic states of the reactants and products are specified.

Figure 6: The TEST Specific, Closed Process, Premixed, Reaction daemon using the IG Mixture Model. The figure shows the state of a reactant mixture of a stoichiometric amount of methane and oxygen gas at standard conditions (1 bar, 298 K).

After state computation, one can then perform a process analysis by either specifying or computing heat transfer, boundary work, boundary temperature, entropy generation, or the change in internal energy or entropy resulting from a closed process from one state to another where a change in energy occurs through heat transfer, work transfer, or both. For example, Figure 6 shows the state of a reactant mixture of a stoichiometric amount of methane and oxygen gas at standard conditions of 1 bar and 298 K. Figure 7 shows the corresponding state of the product mixture of water vapor and carbon dioxide gas. The computed adiabatic flame temperature is shown in the cyan box labeled T2.

Figure 7: The state of the product mixture of water vapor and carbon dioxide gas. The computed adiabatic flame temperature is shown in the cyan box labeled T2.

The IG and IGE models in TEST can be used to perform fairly advanced parametric studies. To illustrate this capability, let us extend the previous example a little further to investigate the effect of oxygen level on the adiabatic flame temperature of methane combustion in air using both the TEST IGE daemon and the Combustion Chamber Simulator RIA [30]. The complete reaction for the IG model is The IGE daemon can be launched by navigating to the TEST Map page and selecting the Closed, Unsteady, Specific, Combustion and Equilibrium Daemon page and then selecting the Chemical Equilibrium (IGE) Model daemon. In the IGE daemon, one defines an initial reactant mixture composition (by mass, volume, or moles) and the allowable product species (but not the product species stoichiometry like is needed in the premixed combustion daemon). For a given pressure and either temperature or enthalpy, the IGE daemon calculates the equilibrium composition and the complete products state, including the equilibrium flame temperature if enthalpy is specified, by minimizing the mixture’s Gibbs function. For example, to define reaction (61) in the IGE daemon, configure the reactant mixture as State 1 by choosing one kmole of CH₄, two kmoles of O₂, and 7.52 kmoles of N₂ in the Composition Panel. Then navigate to the State Panel and set and . Click the Calculate button to find the specific enthalpy of the reactants as shown in Figure 8.

Figure 8: To compute the adiabatic flame temperature of reaction (61) using TEST’s IGE daemon, first define State 1 as the reactant mixture using the Composition Panel and then click the Calculate button in the State Panel to compute the reactant’s specific enthalpy .

Once the specific enthalpy of the reactants is known, the IG model can be used to find the adiabatic flame temperature by defining State 2 to be the products mixture assuming complete combustion. In the Composition Panel, select only the products given in reaction (61) and then return to the State Panel and click the greater than symbol to the right of the state pull down menu to define State 2. Then set and to constrain pressure and specific enthalpy between the two states. Click the Calculate button again and the computed temperature displayed in the cyan box will be the adiabatic flame temperature computed using the IG model as shown in Figure 9.

Figure 9: The IG model can be used within the TEST IGE daemon by choosing only those products specified in the complete reaction (61). The IG model adiabatic flame temperature is computed and displayed in the T2 cyan box.

Finally, to use the IGE model in the IGE daemon, return to the Composition Panel, and pull down the menu labeled Choose Preconfigured Products and select the option Select Common Combustion Species. Then return to the State Panel and click the Calculate button once more. The computed temperature displayed in the cyan box will be the adiabatic flame temperature computed using the IGE model as shown in Figure 10. Notice the IGE model computes a lower adiabatic flame temperature of 2223 K compared with 2326 K as found by the IG model.

Figure 10: The IGE model can be used within the TEST IGE daemon by selecting common combustion species using the preconfigured products pull-down menu in the Composition Panel. The IGE model computes a lower adiabatic flame temperature than the IG model.

All of the screen images of daemons shown in Figures 1 and 6 through Figure 10 are implemented as Java applets. To run any of these daemons, a user must have the Java runtime environment installed on their client PC (available free via URL http://www.java.com/).

Equilibrium distribution for an adiabatic combustion chamber is calculated using the IGE daemon for a reactant mixture of one mole of methane and 2 moles of air. For a stoichiometric amount of oxygen, the IG model computes an adiabatic flame temperature of 2326 K while the IGE model computes the lower value of 2223 K. The IGE result corresponds exactly with published data on flame temperatures [40]. A lower flame temperature found by the IGE model in the higher O₂ mole fraction range is expected since species dissociation will lower flame temperature, as predicted by Le Chatelier’s principle discussed earlier.

The results of a parametric study comparing adiabatic flame temperature versus oxygen reactant mole fraction using the IG and IGE models is shown in Figure 11. The plot clearly shows how dissociation remarkably affects flame temperature.

Figure 11: Adiabatic flame temperature versus oxygen reactant mole fraction. A comparison of IG and IGE model results.

One of the useful features of TEST’s IGE model daemons is the ability to compute the partial pressure equilibrium constant of a reaction at different temperatures. IGE model-based daemons will report ln in the upper right corner of the daemon’s Composition Panel as shown in Figure 12. The IGE model computes the change in the standard Gibbs energy between the reaction products and reactants, . From this value, ln is easily found at a particular temperature by Figure 13 shows a plot comparing ln versus the reciprocal of temperature for the equilibrium reaction The IG and IGE models can also be used to investigate the difference in adiabatic flame temperature versus equivalence ratio which is the ratio of the fuel-to-oxidizer ratio to the stoichiometric fuel-to-oxidizer ratio, where can be mass or number of moles. A value of is at stoichiometry while rich mixtures are greater than 1 and lean mixtures are less than 1. Figure 14 shows a comparison of the IG and IGE models for the constant pressure combustion of propane gas at different equivalence ratios. For the IG model, the complete reaction is For the IGE model we allowed CO, H₂, and OH to form as products which are the predominant product species of combustion other than H₂O and CO₂. The IGE model is shown to compute a lower flame temperature for rich mixtures.

Figure 12: TEST IGE model daemons will compute and report the log of the partial pressure equilibrium constant of a reaction in the upper right corner of the Composition Panel.

Figure 13: Log of the partial pressure equilibrium constant   at specific temperatures for reaction (63). IGE model values agree with published tabulated data [25] to within an error norm of 0.03.

Figure 14: Adiabatic flame temperature versus equivalence ratio for the combustion of propane gas in air at constant pressure.

To assess the performance difference between the iterative and direct methods employed by the IGE model we computed the adiabatic flame temperature of methane and oxygen gas at a standard pressure of 1 bar. The overall reaction is where products varies according the list of species shown in each row of Table 2. As one can see in Table 1, both methods compute the same flame temperature to an accuracy of 1 K. However, the time taken to compute the flame temperature directly is drastically lower than using an iterative approach.

Table 1: Performance comparison of the IGE model iterative and direct methods on the adiabatic flame temperature of methane and oxygen at 1 bar.

Number of Product Species
Iterative Method Wall Time (s)
Direct Method Wall Time (s)
Number of Iterations
Iterative Method Flame Temperature (K)
Direct Method Flame Temperature (K)
Difference

2
281.80
3.19
30
5166
5166
0

3
322.42
15.19
30
5166
5166
0

4
216.17
3.33
21
3842
3842
0

5
92.78
3.70
8
3376
3376
0

6
455.58
4.47
41
3262
3262
0

7
2429.08†
5.05
297
3122
3122
0

8
*
5.67
*
*
3049
*

9
*
7.91
*
*
3049
*

†Run on an Intel 2x Quad Core system (2 Xeon E5530 @ 2.40 GHz CPUs). The other tests were run on an Intel 2x single core system (2 Xeon Nocona @ 2.80 GHz CPUs).

Table 2: Set of possible product species for each adiabatic flame temperature calculation shown in Table 1.

Number of Product Species
Species Set

2
CO₂, H₂O

3
CO₂, H₂O, CO

4
CO₂, H₂O, CO, O₂

5
CO₂, H₂O, CO, O₂, OH

6
CO₂, H₂O, CO, O₂, OH, O

7
CO₂, H₂O, CO, O₂, OH, O, H₂

8
CO₂, H₂O, CO, O₂, OH, O, H₂, H

9
CO₂, H₂O, CO, O₂, OH, O, H₂, H, O₃

Using the IGE model to perform a type computation can be seen in investigating the effect to which chemical equilibrium plays a role in air compressor efficiency. To illustrate, suppose an air compressor compresses incoming air, assumed to be a mixture of 21% O₂ and 79% N₂ by volume, from the standard state of 298 K and 1 bar, to an exit state of 800 kPa and variable temperature . For values of , where is the exit temperature reached if the inlet air was compressed reversibly and adiabatically (i.e., isentropically), we can compute the compressor isentropic efficiency, , using the IG and IGE model within the TEST Open System, Steady-State, Single Flow Equilibrium daemon.

The compressor isentropic efficiency is defined as where is the specific enthalpy of the exit state for an isentropic compression and is the specific enthalpy of the exit state for an actual compression. The term can be computed using the daemon by assigning the inlet and outlet pressures accordingly and setting the specific entropy of both states to be equal. The daemon uses the iterative scheme to compute the exit temperature based on an isentropic compression process. For this particular example problem the isentropic, or ideal, exit temperature to compress air in a standard state to 800 kPa is 535 K. To obtain the for a particular exit temperature using the IG model we simply use the daemon to find the exit state specific enthalpy of a reactant mixture of 1 kmole of O2 and 3.76 kmols of N₂. For the IGE model we use the daemon to find the exit state specific enthalpy of air products which consists of some distribution of N, N₂, N₂O, N₂O₄, NO, NO₂, O, O₂, and O₃ in equilibrium. These nine species represent the feasible species that might form as air is compressed to a high pressure and temperature. Figure 15 shows a plot of versus compressor exit temperature. From the plot one can see that chemical equilibrium does not have an appreciable effect on compressor efficiency, except at very high exit temperatures where the discrepancy is only slight. At low temperatures air will not dissociate and so the IG model matches the IGE model. At very high temperatures >3000 K air will begin to dissociate and form oxides of nitrogen as well as monatomic oxygen and nitrogen gas. At 6000 K air is 31% monatomic oxygen, 6% monatomic nitrogen, and 2% nitric oxide, with the remainder being mostly diatomic nitrogen gas.

Figure 15: Investigating the effect to which chemical equilibrium plays a role in air compressor efficiency. The IG model shows a slight deviation from the IGE model only at very high exit temperatures.

As a final interesting example of how the TEST IGE model can provide a more accurate representation of equilibrium effects, Figure 16 shows a plot of the molar mass of the products of hydrogen combustion with oxygen versus the ratio of oxygen to hydrogen fuel in a rocket combustion chamber. The TEST IG model solves the complete reaction while the IGE model solves the equilibrium reaction with nine product species, From the plot one can see how dissociation, association, and recombination lower the exhaust molecular weight at higher oxygen to hydrogen ratios. Because a rocket must internally carry its oxidizer and fuel, minimizing the total mass of oxidizer and fuel increases forward acceleration. The molecular weight of oxygen is sixteen-times greater than that of hydrogen and so overall rocket motor performance is improved at lower oxidizer/fuel ratios. For a fixed mass of oxidizer and fuel, the change in forward momentum can only be increased by an equal and opposite change in exhaust momentum. The IGE model clearly shows, however, that equilibrium plays a role in reducing exhaust momentum at higher oxidizer/fuel ratios since the molecular weight of the exhaust is less than in an ideal case of complete combustion.

Figure 16: Product mixture molar mass versus the ratio of oxygen to hydrogen fuel in a rocket combustion chamber. The TEST IG model solves the complete reaction (68) while the IGE model solves the equilibrium reaction (69).

5. Conclusion

TEST is a collection of thermodynamic databases, software, and courseware, accessible on the Internet from http://www.thermofluids.net/. TEST is freely accessible from any academic institution and is currently used in more than 100 institutions around the world. The thermodynamic calculators in TEST, called daemons, have been developed using an object-oriented paradigm in Java so that code can be reused. All TEST daemons have been thoroughly tested under different versions of browsers such as IE, Firefox, Chrome, and Safari running on Windows, MacOS, and Linux platforms. The only browser plug-in required to run TEST daemons is version 4 (or greater) of the Java Platform.

Gas mixtures are modeled in TEST by the perfect gas or ideal gas mixture models. In this work we extended the ideal gas (IG model) mixture model into an ideal gas equilibrium (IGE model) mixture model by incorporating chemical equilibrium calculations as part of the state evaluation. Using TEST’s intuitive graphical interface, users can define a reaction by specifying the composition of a reacting mixture and the possible species that may form as products, and compute a thermodynamic state that is based on an equilibrium distribution. For a given pressure and either temperature, enthalpy, or entropy, a product mixture’s Gibbs function is minimized subject to atomic mass and possibly either energetic or entropic constraints. The resulting equilibrium composition along with thermodynamic properties of the mixture are calculated and displayed.

Whereas the IG model computes state variables based on a frozen mixture, the IGE model yields a more accurate determination of mixture properties such as specific enthalpy, specific entropy, and temperature. Thus the IGE model is well suited for computing the equilibrium constant of a reaction or final temperatures resulting from adiabatic or isentropic processes. The differences between the IG and IGE models become apparent at high temperatures where species dissociation occurs according to Le Chatelier’s Principle. The implementation of the IGE model in TEST uses two approaches when solving and problems: a direct method and an iterative method. Both methods are accurate in that they compute equivalent states. The iterative method is simpler to comprehend but is not practical to use when computing the equilibrium distribution for mixtures containing many species as the number of required iterations increases greatly and consequently the overall time to achieve convergence becomes unacceptable. The approach taken by TEST is unique in that equilibrium computations are performed in the background, without requiring any major change in the familiar interface used in other state daemons. Thus where the IG model can be used, the IGE can be used with ease to obtain a more accurate state. Thermodynamic states and distributions calculated by the TEST IGE model agree well with results from other established chemical equilibrium applications such as NASA CEA [11, 36], and STANJAN [39].

Nomenclature

: | Number of different elements in the reacting system
: | Number of atoms of type in ideal gas
: | Number of atoms of type in the reacting system
: | Standard state specific molar heat capacity of species , J mol⁻¹ K⁻¹
: | Standard state heat of formation
: | Index used to reference a gaseous species
: | Partial pressure equilibrium constant
: | Index used to reference a condensed species
: | Number of unique species in the product composition
: | Number of moles of species
: | Total number of moles in the equilibrium composition
: | Species in the gas phase
: | Species in the liquid phase
: | Enthalpy, kJ
: | Specific molar enthalpy of the reactant mixture, J mol⁻¹
: | Standard state specific molar enthalpy of species , J mol⁻¹
: | Temperature, K
: | Entropy, KJ·K⁻¹
: | Specific molar entropy of the reactant mixture, J mol−1 K⁻¹
: | Specific molar enthalpy of species , J mol−1 K⁻¹
: | Standard state specific molar entropy of species , J mol−1 K⁻¹
: | Internal energy, kJ
: | Pressure, atm
: | Partial pressure of species , atm
: | Standard state pressure = 1 atm
: | Ideal gas constant = 8.314472 J·K⁻¹·mol⁻¹
: | Volume,
: | Mole fraction of species
: | Compressor isentropic efficiency
: | Equivalence ratio
: | Chemical potential of the th species, J·mol⁻¹
: | Standard state chemical potential of the th species, J·mol⁻¹.

Acknowledgment

The authors gratefully acknowledge support from NSF through the Cyberinfrastructure CI-TEAM Grant 0753283 and support from NASA.

References

P. Atkins and J. de Paula, Physical Chemistry, W.H. Freeman, New York, NY, USA, 8th edition, 2006.

O. Redlich and J. N. S. Kwong, “On the thermodynamics of solutions,” Chemical Reviews, vol. 44, no. 1, pp. 233–244, 1949.

D. Y. Peng and D. B. Robinson, “A new two-constant equation of state,” Industrial and Engineering Chemistry Fundamentals, vol. 15, no. 1, pp. 59–64, 1976.

D. M. Kerrick and G. K. Jacobs, “A modified Redlich-Kwong equation for H₂O, CO₂, and H₂O-CO₂ mixtures at elevated pressures and temperatures,” American Journal of Science, vol. 281, no. 6, pp. 735–767, 1981.

W. B. White, S. M. Johnson, and G. B. Dantzig, “Chemical equilibrium in complex mixtures,” The Journal of Chemical Physics, vol. 28, no. 5, pp. 751–755, 1958.

F. J. Zeleznik and S. Gordon, “An analytical investigation of three general methods of calculating chemical-equilibrium compositions,” Tech. Rep. NASA TN D-473, 1960.

F. J. Zeleznik and S. Gordon, “Simultaneous least-squares approximation of a function and its first integrals with application to thermodynamic data,” Tech. Rep. NASA TN D-767, NASA, New York, NY, USA, 1961.

F. J. Zeleznik and S. Gordon, “A general IBM 704 or 7090 computer program for computation of chemical equilibrium compositions, rocket performance, and Chapman-Jouguet Detonations,” Tech. Rep. NASA TN D-1454, NASA, New York, NY, USA, 1962.

F. J. Zeleznik and S. Gordon, “Calculation of detonation properties and effect of independent parameters on gaseous detonations,” American Rocket Society Journal, vol. 32, no. 4, pp. 606–615, 1962.

F. J. Zeleznik and S. Gordon, “Calculation of complex chemical equilibria,” Industrial & Engineering Chemistry, vol. 60, no. 6, pp. 27–57, 1968.

S. Gordon and B. J. McBride, “Computer program for calculation of complex chemical equilibrium compositions and applications—I. analysis,” Tech. Rep. 1311, NASA Glenn Research Center, Cleveland, Ohio, USA, 1994.

M. Castier, P. Rasmussen, and A. Fredenslund, “Calculation of simultaneous chemical and phase equilibria in nonideal systems,” Chemical Engineering Science, vol. 44, no. 2, pp. 237–248, 1989.

M. L. Michelsen, “The isothermal flash problem. Part I. Stability,” Fluid Phase Equilibria, vol. 9, no. 1, pp. 1–19, 1982.

M. L. Michelsen, “The isothermal flash problem. Part II. Phase-split calculation,” Fluid Phase Equilibria, vol. 9, no. 1, pp. 21–40, 1982.

I. G. Economou, M. D. Donohue, and L. W. Hunter, “Equation-of-state calculations of chemical reaction equilibrium in nonideal systems,” International Journal of Thermophysics, vol. 14, no. 2, pp. 199–213, 1993.

G. Soave, “Equilibrium constants from a modified Redlich-Kwong equation of state,” Chemical Engineering Science, vol. 27, no. 6, pp. 1197–1203, 1972.

R. Stryjek and J. H. Vera, “PRSV2: a cubic equation of state for accurate vapor—liquid equilibria calculations,” The Canadian Journal of Chemical Engineering, vol. 64, no. 5, pp. 820–826, 1986.

P. Proust and J. H. Vera, “PRSV: the stryjek-vera modification of the peng-robinson equation of state. Parameters for other pure compounds of industrial interest,” Canadian Journal of Chemical Engineerin, vol. 67, no. 1, pp. 170–173, 1989.

R. Stryjek and J. H. Vera, “PRSV: an improved peng—robinson equation of state for pure compounds and mixtures,” Canadian Journal of Chemical Engineering, vol. 64, no. 2, pp. 323–333, 1986.

R. Stryjek and J. H. Vera, “PRSV—an improved peng-Robinson equation of state with new mixing rules for strongly nonideal mixtures,” Canadian Journal of Chemical Engineering, vol. 64, no. 2, pp. 334–340, 1986.

M. Llano-Restrepo and Y. M. Muñoz-Muñoz, “Combined chemical and phase equilibrium for the hydration of ethylene to ethanol calculated by means of the Peng-Robinson-Stryjek-Vera equation of state and the Wong-Sandler mixing rules,” Fluid Phase Equilibria, vol. 307, no. 1, pp. 45–57, 2011.

D. S. Abrams and J. M. Prausnitz, “Statistical thermodynamics of liquid mixtures: a new expression for the excess Gibbs energy of partly or completely miscible systems,” AIChE Journal, vol. 21, no. 1, pp. 116–128, 1975.

S. I. Sandler, Chemical, Biochemical, and Engineering Thermodynamics, John Wiley & Sons, New York, NY, USA, 4th edition, 2006.

J. A. Beattie and O. C. Bridgeman, “A new equation of state for fluids. I. Application to gaseous ethyl ether and carbon dioxide,” Journal of the American Chemical Society, vol. 49, no. 7, pp. 1665–1667, 1927.

Y. A. Cengel and M. A. Boles, Thermodynamics—An Engineering Approach, McGraw–Hill, New York, NY, USA, 4th edition, 2002.

M. Benedict, G. B. Webb, and L. C. Rubin, “An empirical equation for thermodynamic properties of light hydrocarbons and their mixtures: I. Methane, ethane, propane and n-butane,” The Journal of Chemical Physics, vol. 8, no. 4, pp. 334–345, 1940.

H. W. Cooper and J. C. Goldfrank, “B-W-R constants and new correlations,” Hydrocarbon Processing, vol. 46, no. 12, pp. 141–146, 1967.

S. Bhattacharjee, The Expert System for Thermodynamics: A Visual Tour, Prentice Hall, New York, NY, USA, 2002.

S. Bhattacharjee, “TEST—the expert system for thermodynamics,” in Proceedings of the American Society for Engineering Education Annual Conference & Exposition, p. 15, Nashville, Tenn, USA, 2003.

M. Patterson, C. Paolini, and S. Bhattacharjee, “Design and implementation of a rich internet application (RIA) for the simulation of a combustion chamber,” in Proceedings of the ASEE Annual Conference and Exposition, pp. 1–16, Louisville, Ky, USA, June 2010.

S. Bhattacharjee and C. P. Paolini, “Property evaluation in The Expert System for Thermodynamics ("TEST") web application,” Journal of Computer Coupling of Phase Diagrams and Thermochemistry, vol. 33, no. 2, pp. 343–352, 2009.

S. Bhattacharjee and C. P. Paolini, “The chemical thermodynamic module of the expert system for thermodynamics (“TEST”) web application,” in Proceedings of the ASEE Annual Conference and Exposition, pp. 1–17, Austin, Tex, USA, June 2009.

C. Paolini, K. Bobba, P. Surana, and S. Bhattacharjee, “A Java based web application for performing chemical equilibrium analysis in thermodynamics courses,” in Proceedings of the 36th ASEE/IEEE Frontiers in Education Conference (FIE '06), San Diego, Calif, USA, October 2006.

C. P. Paolini and S. Bhattacharjee, “An object-oriented online tool for solving generalized chemical equilibrium problems,” in Proceedings of the ASME International Mechanical Engineering Congress and Exposition (IMECE '08), Boston, Mass, USA, October/November 2008.

C. Paolini and S. Bhattacharjee, “Solving chemical equilibrium problems online,” Journal of Chemical Education, vol. 87, no. 4, p. 456, 2010.

B. J. McBride and S. Gordon, “Computer program for calculation of complex chemical equilibrium compositions and applications—II. users manual and program description,” Tech. Rep. 1311-P2, NASA Glenn Research Center, Cleveland, Ohio, USA, 1996.

S. Gordon and B. J. McBride, “Computer program for computation of complex chemical equilibrium compositions, rocket performance, incident and reflected shocks, and Chapman-Jouguet detonations,” Tech. Rep. SP-273, NASA, New York, NY, USA, 1971.

S. Gordon, F. J. Zeleznik, and V. N. Huff, “A general method for automatic computation of equilibrium compositions and theoretical rocket performance of propellants,” Tech. Rep. TN D-132, NASA, New York, NY, USA, 1959.

W. C. Reynolds, “The element potential method for chemical equilibrium analysis: implementation in the interactive program STANJAN version 3,” Tech. Rep., Stanford University, Department of Mechanical Engineering, Stanford, Calif, USA, 1986.

V. Babrauskas, Temperatures in Flames and Fires, Fire Science and Technology, Washington, DC, USA, 2006.