About this Journal Submit a Manuscript Table of Contents
Journal of Thermodynamics
Volume 2012 (2012), Article ID 870631, 18 pages
http://dx.doi.org/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

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, NH3, 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,𝑝𝑣=𝑅𝑇.(1) 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𝑝𝑣=𝑧𝑅𝑇,(2) 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𝑃=𝑅𝑇𝑣𝑏𝑎𝑣2.(3) 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𝜕𝑃𝜕𝑣𝑇𝑐=0,𝜕2𝑃𝜕𝑣2𝑇𝑐=0,(4) 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],𝑃=𝑅𝑇𝑣𝑏𝑎𝑇𝑣(𝑣+𝑏).(5) 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 [610]. 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𝑃=𝑅𝑇𝑣𝑏𝑎𝛼𝑣(𝑣+𝑏),(6) 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𝑃=𝑅𝑇𝑣𝑏𝑎𝛼𝑣(𝑣+𝑏)+𝑏(𝑣𝑏),(7) where three parameters, 𝜅1, 𝜅2, and 𝜅3 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 𝜅1, 𝜅2, and 𝜅3 for many fluids are found in the literature [1820]. 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𝑃=𝑅𝑢𝑇𝑣21𝑐𝑣𝑇3𝑣+𝑏𝐴𝑣2,𝐴=𝐴01𝑎𝑣,𝐵=𝐵01𝑏𝑣,(8) where 𝐴 and 𝐵 are empirical parameters, it is fairly accurate for species with densities up to 0.8𝜌cr [25]. The Benedict-Webb-Rubin EOS [26],𝑃=𝑅𝑢𝑇𝑣+𝐵0𝑅𝑢𝑇𝐴0𝐶0𝑇21𝑣2+𝑏𝑅𝑢𝑇𝑎𝑣3+𝑎𝛼𝑣6+𝑐𝑣3𝑇21𝛾𝑣2𝑒𝛾/𝑣2,(9) where the constants 𝑎, 𝐴0, 𝑏, 𝐵0, 𝑐, 𝐶0, 𝛼, and 𝛾 have been tabulated for many gas species [27], it is accurate for species with densities up to 2.5𝜌cr. 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 [3234] 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., 𝑐𝑝=𝑐𝑣=constant). 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 𝑝1=1atm and 𝑇1=3000K. The notation 𝑝1 and 𝑇1 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 O2 molecules, and O atoms will form and this mixture will have different thermodynamic properties than a pure amount of O2 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.

870631.fig.001
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.
870631.fig.002
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 O2 molecules and O atoms are plotted with respect to temperature. From this figure one can see that at 6000 K a sample of O2 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 O2molecules 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 O2 molecule will dissociate into two O atoms. The bond dissociation energy D°  for an O=O bond which is broken through the reaction O2O+O is defined as the standard state enthalpy change for the reaction at a specified temperature,𝐷=2Δ𝐻𝑓(O)Δ𝐻𝑓O2.(10) 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 O2 and O for a 1 mole initial mixture of O2 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 O2 bonds are broken to form oxygen atoms is matched by the rate at which two oxygen atoms recombine to form an oxygen molecule.O22O.(11) The amount of O2 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 O2 molecules that dissociate into 2𝛼 moles of O atoms, leaving 1-𝛼 moles of O2, we can express 𝐾𝑝 for reaction (11) in terms of the partial pressures of the reactant O2 and the product O,𝐾𝑝=𝑝2O𝑝O2=𝑦2O(𝑃)2𝑦O2𝑃=[2𝛼/((1𝛼)+(2𝛼))]2(1𝛼)/((1𝛼)+(2𝛼))𝑃=4𝛼21𝛼2𝑃.(12) 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 O2molecules. As the number of O atoms increases, the partial pressure of O atoms increases and this increases the equilibrium constant 𝐾𝑝.

870631.fig.003
Figure 3: The mole fractions of O2 molecules and O atoms in a one mole initial sample of oxygen gas at standard pressure. As temperature increases, O2 molecules dissociate into O atoms. At 6000 K, a one mole initial sample of oxygen gas will contain 2 moles of O atoms.
870631.fig.004
Figure 4: Log of the partial pressure equilibrium constant 𝐾𝑝 versus the reciprocal of temperature for reaction (11).

Using the van’t Hoff equation,𝜕ln𝐾𝑝𝜕(1/𝑇)𝑝=Δ𝐻𝑅,(13) 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 O2 molecules and O atoms are present.

870631.fig.005
Figure 5: Standard state molar enthalpy of reaction Δ𝐻 versus temperature for reaction (11). Δ𝐻>0 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𝐺=𝐻𝑇𝑆=(𝑈+𝑃𝑉)𝑇𝑆=𝑚𝑗=1𝜇𝑗𝑛𝑗.(14) A system will be in equilibrium when the change in Gibbs free energy of the system vanishes to zero.

Differentiating (14) we obtain𝑑𝐺=𝑑𝑈+𝑃𝑑𝑉+𝑉𝑑𝑃𝑇𝑑𝑆𝑆𝑑𝑇.(15) At a constant temperature 𝑇 and pressure 𝑃, (15) reduces to𝑑𝐺=𝑑𝑈+𝑃𝑑𝑉𝑇𝑑𝑆.(16) The combined first and second laws of thermodynamics requires𝑑𝑈𝑇𝑑𝑆+𝑃𝑑𝑉0𝑑𝐺0,(17) 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 𝑗,𝜇𝑗=𝜇𝑗+𝑅𝑇ln𝑃𝑗𝑃.(18) Substituting (18) into (14) and dividing by 𝑅𝑇 to obtain a dimensionless form, we obtain𝐺𝑅𝑇=𝑚𝑗=1𝑛𝑗𝜇𝑗𝑅𝑇+𝑛𝑗ln𝑃𝑗𝑃,(19) butln𝑃𝑗𝑃=ln𝑃𝑗𝑃𝑃𝑃=ln𝑦𝑗𝑃𝑃𝑃𝑃=ln𝑦𝑗+ln𝑃𝑃,(20) where the last term, ln(𝑃/𝑃), reduces to a constant since the system total pressure 𝑃 is given as an input to the problem. Using (20), (19) now becomes𝐺𝑅𝑇=𝑚𝑗=1𝑛𝑗𝜇𝑗𝑅𝑇+𝑛𝑗ln𝑦𝑗+𝑛𝑗ln𝑃𝑃.(21) The minimum stationary point of (14) will be the vector of species molar values 𝑛where 𝑑𝐺vanishes. Differentiating 𝐺 in (14) we obtain𝑑𝐺=𝑚𝑗=1𝑛𝑗𝑑𝜇𝑗+𝑚𝑗=1𝜇𝑗𝑑𝑛𝑗,(22) but from the isothermal, isobaric Gibbs-Duhem equation we know that𝑚𝑗=1𝑛𝑗𝑑𝜇𝑗=0,(23) and so we seek the unique vector 𝑛such that𝑚𝑗=1𝜇𝑗𝑅𝑇+ln𝑛𝑗ln𝑛+ln𝑃𝑃𝑑𝑛𝑗=0.(24) 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𝑚𝑗=1𝐴𝑖,𝑗𝑛𝑗=𝑏𝑖,𝑖=1,,𝑎.(25) 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𝑛𝑗0,𝑗,1𝑗𝑚.(26) We can express the atomic constraint equations given by (25) and (26) as equality constraints𝜙𝑖=𝜙𝑖𝑛1,𝑛2,,𝑛𝑚=𝑚𝑗=1𝐴𝑖,𝑗𝑛𝑗𝑏𝑖=0.(27) There will be 𝑎 equations of the form (27) for each atom 𝑖 present in the system. We define the Lagrangian function =(𝑛,𝜆) as 𝑛,𝜆=𝜇𝑛+𝑎𝑖=1𝜆𝑖𝜙𝑖.(28) The solution of (24) is thus the minimum stationary point where the gradient of (28) vanishes, that is,𝑛,𝜆=𝑚𝑗=1𝜕𝜇𝜕𝑛𝑗𝑝,𝑇,𝑛𝑖𝑗𝑑𝑛𝑗+𝑎𝑖=1𝜆𝑖𝑚𝑗=1𝜕𝜙𝑖𝜕𝑛𝑗𝑑𝑛𝑗=0.(29) The TEST IGE model uses a Newton-Raphson technique to solve (29). By expanding (29) we can write𝑛,𝜆=𝑚𝑗=1𝜕𝜇𝜕𝑛𝑗𝑝,𝑇,𝑛𝑖𝑗𝑑𝑛𝑗+𝑎𝑖=1𝜆𝑖𝑚𝑗=1𝜕𝜙𝑖𝜕𝑛𝑗𝑑𝑛𝑗=𝑚𝑗=1𝜇𝑗𝑑𝑛𝑗+𝑎𝑖=1𝜆𝑖𝜕𝜙𝑖𝜕𝑛1𝑑𝑛1+𝜕𝜙𝑖𝜕𝑛2𝑑𝑛2++𝜕𝜙𝑖𝜕𝑛𝑚𝑑𝑛𝑚=𝑚𝑗=1𝜇𝑗𝑑𝑛𝑗+𝜆1𝜕𝜙1𝜕𝑛1𝑑𝑛1+𝜕𝜙1𝜕𝑛2𝑑𝑛2++𝜕𝜙1𝜕𝑛𝑚𝑑𝑛𝑚+𝜆2𝜕𝜙2𝜕𝑛1𝑑𝑛1+𝜕𝜙2𝜕𝑛2𝑑𝑛2++𝜕𝜙2𝜕𝑛𝑚𝑑𝑛𝑚++𝜆𝑎𝜕𝜙𝑎𝜕𝑛1𝑑𝑛1+𝜕𝜙𝑎𝜕𝑛2𝑑𝑛2++𝜕𝜙𝑎𝜕𝑛𝑚𝑑𝑛𝑚=𝜇1+𝜆1𝜕𝜙1𝜕𝑛1+𝜆2𝜕𝜙2𝜕𝑛1++𝜆𝑎𝜕𝜙𝑎𝜕𝑛1𝑑𝑛1+𝜇2+𝜆1𝜕𝜙1𝜕𝑛2+𝜆2𝜕𝜙2𝜕𝑛2++𝜆𝑎𝜕𝜙𝑎𝜕𝑛2𝑑𝑛2++𝜇𝑚+𝜆1𝜕𝜙1𝜕𝑛𝑚+𝜆2𝜕𝜙2𝜕𝑛𝑚++𝜆𝑎𝜕𝜙𝑎𝜕𝑛𝑚𝑑𝑛𝑚=0.(30) 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𝜇𝑗+𝑎𝑖=1𝜆𝑖𝜕𝜙𝑖𝜕𝑛𝑗=0,(31) for all 𝑗. By observing𝜆𝑖𝜕𝜙𝑖𝜕𝑛𝑗=𝜆𝑖𝐴𝑖,𝑗,(32) solving (30) amounts to solving a system of 𝑚 equations for each species 𝑗 where the 𝑗th species equation is given by𝜇𝑗+𝑎𝑖=1𝜆𝑖𝐴𝑖,𝑗=0,(33) and additionally, 𝑎 population constraint equations for each element 𝑖 where the 𝑖th constraint equation is given by𝑚𝑗=1𝐴𝑖,𝑗𝑛𝑗𝑏𝑖=0,(34) and finally one constraint equation governing the total number of moles present in the system,𝑚𝑗=1𝑛𝑗𝑛=0.(35) From (33), (34), and (35), we must solve a system of 𝑚+𝑎+1equations in 𝑎+𝑚+1unknowns. The IGE model uses the Newton-Raphson method to iteratively solve this system of equations. To construct our system matrix we let 𝑓𝑗=𝑓𝑗(𝑥)=𝑓𝑖(𝑥1,𝑥2,,𝑥𝑛) represent the 𝑗th species equation of (33), or𝑓𝑗=𝜇𝑗+𝑎𝑖=1𝜆𝑖𝐴𝑖,𝑗=0.(36) Expanding 𝜇𝑗 using (21) we have𝑓𝑗=𝜇𝑗𝑅𝑇+ln𝑛𝑗ln𝑛+ln𝑃𝑃+𝑎𝑖=1𝜆𝑖𝑅𝑇𝐴𝑖,𝑗=0.(37) By letting𝑥1=ln𝑛𝑗,𝑥2=ln𝑛,𝑥3=𝜆1𝑅𝑇,𝑥4=𝜆2𝑅𝑇,,𝑥𝑎+2=𝜆𝑎𝑅𝑇,(38) we will obtain𝜕𝑓𝑗𝜕𝑥1=𝜕𝑓𝑗𝜕ln𝑛𝑗=1,𝜕𝑓𝑗𝜕𝑥2=𝜕𝑓𝑗𝜕(ln𝑛)=1,𝜕𝑓𝑗𝜕𝑥3=𝜕𝑓𝑗𝜕𝜆1/𝑅𝑇=𝐴1,𝑗,,(39) and our Newton-Raphson equation becomes𝑛𝑖=1𝜕𝑓𝜕𝑥𝑖𝛿𝑥𝑖=𝑓𝑥𝛿ln𝑛𝑗𝛿ln𝑛+𝑎𝑖=1𝐴𝑖,𝑗𝜆𝑖𝑅𝑇=𝜇𝑗𝑅𝑇𝑎𝑖=1𝜆𝑖𝑅𝑇𝐴𝑖,𝑗.(40) We will have 𝑚 equations of type (40) for each species 𝑗. For the population constraint equations given by (34), we let 𝑓𝑖=𝑓𝑖(𝑥)=𝑓𝑖(𝑥1,𝑥2,,𝑥𝑛) represent the 𝑖th population constraint equation such that𝑓𝑖=𝑚𝑗=1𝐴𝑖,𝑗𝑛𝑗𝑏𝑖=𝑚𝑗=1𝐴𝑖,𝑗expln𝑛𝑗𝑏𝑖=0.(41) Then𝜕𝑓𝑖𝜕𝑥𝑗=𝜕𝑓𝑖𝜕ln𝑛𝑗=𝐴𝑖,𝑗𝑛𝑗,(42) and our Newton equation for (34) becomes𝑚𝑗=1𝐴𝑖,𝑗𝑛𝑗𝛿ln𝑛𝑗=𝑏𝑖𝑚𝑗=1𝐴𝑖,𝑗𝑛𝑗.(43) 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 𝑓=𝑓(𝑥)=𝑓(𝑥1,𝑥2,,𝑥𝑛) represent (35),𝑓=𝑚𝑗=1expln𝑛𝑗exp(ln𝑛)=0,(44) where𝜕𝑓𝜕𝑥𝑗=𝜕𝑓𝜕ln𝑛𝑗=𝑛𝑗,𝜕𝑓𝜕(ln𝑛)=𝑛,(45) and thus the Newton equation for the total system moles constraint becomes𝑛𝑖=1𝜕𝑓𝜕𝑥𝑖𝛿𝑥𝑖=𝑓𝑥𝑚𝑗=1𝑛𝑗𝛿ln𝑛𝑗𝑛𝛿ln𝑛=exp(ln𝑛)𝑚𝑗=1expln𝑛𝑗=𝑛𝑚𝑗=1𝑛𝑗.(46) We can represent our system of Newton equations in matrix form by combining (40), (43), and (46) to obtain100𝐴1,1𝐴2,1𝐴𝑎,11010𝐴1,2𝐴2,2𝐴𝑎,21001𝐴1,𝑚𝐴2,𝑚𝐴𝑎,𝑚1𝐴1,1𝑛1𝐴1,2𝑛2𝐴1,𝑚𝑛𝑚0000𝐴2,1𝑛1𝐴2,2𝑛2𝐴2,𝑚𝑛𝑚0000𝐴𝑎,1𝑛1𝐴𝑎,2𝑛2𝐴𝑎,𝑚𝑛𝑚0000𝑛1𝑛2𝑛𝑚000𝑛𝛿ln𝑛1𝛿ln𝑛2𝛿ln𝑛𝑚𝜆1𝑅𝑇𝜆2𝑅𝑇𝜆a𝑅𝑇𝛿ln𝑛=𝜇1𝑅𝑇𝑎𝑖=1𝜆𝑖𝐴𝑖,1𝑅𝑇𝜇2𝑅𝑇𝑎𝑖=1𝜆𝑖𝐴𝑖,2𝑅𝑇𝜇𝑚𝑅𝑇𝑎𝑖=1𝜆𝑖𝐴𝑖,𝑚𝑅𝑇𝑏1𝑚𝑗=1𝐴1,𝑗𝑛𝑗𝑏2𝑚𝑗=1𝐴2,𝑗𝑛𝑗𝑏𝑎𝑚𝑗=1𝐴𝑎,𝑗𝑛𝑗𝑛𝑚𝑗=1𝑛𝑗.(47)

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=.(48) The specific enthalpy of the product mixture is determined after each iteration of (47) through the summation=𝑚𝑗=1𝑛𝑗𝑗,(49) 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𝑚𝑗=1𝑛𝑗𝑗𝑅𝑇Δln𝑛𝑗.(50) For an 𝑝 problem, temperature 𝑇varies and approaches the adiabatic flame temperature as (47) converges. As 𝑇 varies, the product mixture enthalpy also varies as𝑚𝑗=1𝑛𝑗𝑐𝑝,𝑗𝑅Δln𝑇,(51) where 𝑐𝑝,𝑗 is the standard state specific molar heat capacity of species 𝑗. Equation (51) is clear sinceΔ=𝑐𝑝Δ𝑇for 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𝑚𝑗=1𝑛𝑗𝑗𝑅𝑇Δln𝑛𝑗Enthalpychangethroughachangeinspeciesmoles+𝑚𝑗=1𝑛𝑗𝑐𝑝,𝑗𝑅Δln𝑇Enthalpychangethroughachangeintemperature+𝑅𝑇Enthalpyoftheproductmixtureatthecurrentiteration=𝑅𝑇Enthalpyofthereactants,aconstant(52) or𝑚𝑗=1𝑛𝑗𝑗𝑅𝑇Δln𝑛𝑗+𝑚𝑗=1𝑛𝑗𝑐𝑝,𝑗𝑅Δln𝑇=𝑅𝑇.(53) Equation (53) is inserted into (47) and repeatedly solved yielding correction terms Δln𝑛𝑗 andΔln𝑇. 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𝑠=𝑠,(54) 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𝑠=𝑚𝑗=1𝑛𝑗𝑠𝑗,(55) where 𝑠𝑗 is the specific molar enthalpy of species𝑗 given by𝑠𝑗=𝑠𝑗𝑅ln𝑛𝑗𝑛𝑅ln𝑃,(56) and 𝑠𝑗 is the standard state specific molar entropy of species 𝑗 given by𝑠=𝑇0𝑐𝑝(𝑇)𝑇𝑑𝑇.(57) As the number of moles of each species 𝑛𝑗 varies, the current value of the product mixture entropy varies as𝑚𝑗=1𝑛𝑗𝑠𝑗𝑅Δln𝑛𝑗.(58) For an 𝑠𝑝 problem, temperature 𝑇 varies and approaches an isentropic temperature as (47) converges. As 𝑇 varies, the product mixture entropy also varies as𝑚𝑗=1𝑛𝑗𝑐𝑝,𝑗𝑅Δln𝑇=𝑚𝑗=1𝑛𝑗Δln𝑃𝑗=𝑚𝑗=1𝑛𝑗Δln𝑛𝑗𝑛+𝑚𝑗=1𝑛𝑗Δln𝑃𝑃=0=𝑚𝑗=1𝑛𝑗Δln𝑛𝑗𝑚𝑗=1𝑛𝑗Δln𝑛=𝑛𝑚𝑗=1𝑛𝑗.(59) The last equality in (59) is from (46). Our additional constraint equation for an isentropic, isobaric process is thus𝑚𝑗=1𝑛𝑗𝑠𝑗𝑅Δln𝑛𝑗+𝑚𝑗=1𝑛𝑗𝑐𝑝,𝑗𝑅Δln𝑇=𝑠𝑠𝑅+𝑛𝑚𝑗=1𝑛𝑗.(60) 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.

alg1
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.

870631.fig.006
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.

870631.fig.007
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 isCH4+2O2+3.76N2CO2+2H2O+7.52N2.(61) 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 CH4, two kmoles of O2, and 7.52 kmoles of N2 in the Composition Panel. Then navigate to the State Panel and set𝑝1=100kPa and 𝑇1=298.15K. Click the Calculate button to find the specific enthalpy of the reactants as shown in Figure 8.

870631.fig.008
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 1.

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 𝑝1=𝑝2and 2=1 to constrain pressure and specific enthalpy between the two states. Click the Calculate button again and the computed temperature displayed in the 𝑇2 cyan box will be the adiabatic flame temperature computed using the IG model as shown in Figure 9.

870631.fig.009
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 𝑇2 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.

870631.fig.0010
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 O2 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.

870631.fig.0011
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 𝑇 byln𝐾𝑝(𝑇)=Δ𝑟𝑥𝑛𝐺𝑅𝑇.(62) Figure 13 shows a plot comparing ln 𝐾𝑝 versus the reciprocal of temperature for the equilibrium reactionCO2CO+12O2.(63) 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,𝜙=𝑚fuel/𝑚oxidizer𝑚fuel/𝑚oxidizerstoichiometric,(64) where 𝑚 can be mass or number of moles. A value of 𝜑=1is 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 isC3H8+5O2+3.76N23CO2+4H2O+18.80N2.(65) For the IGE model we allowed CO, H2, and OH to form as products which are the predominant product species of combustion other than H2O and CO2. The IGE model is shown to compute a lower flame temperature for rich mixtures.

870631.fig.0012
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.
870631.fig.0013
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.
870631.fig.0014
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 isCH4+2O2𝑝𝑟𝑜𝑑𝑢𝑐𝑡𝑠,(66) 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.

tab1
Table 1: Performance comparison of the IGE model iterative and direct methods on the adiabatic flame temperature of methane and oxygen at 1 bar.
tab2
Table 2: Set of possible product species for each adiabatic flame temperature calculation shown in Table 1.

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% O2 and 79% N2 by volume, from the standard state of 298 K and 1 bar, to an exit state of 800 kPa and variable temperature 𝑇2. For values of 𝑇2>𝑇𝑠, 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𝜂𝐶=IsentropiccompressorworkActualcompressorwork=𝑤𝑠𝑤𝑎=Δ𝑠Δ𝑎=2𝑠12𝑎1,(67) where 2𝑠 is the specific enthalpy of the exit state for an isentropic compression and 2𝑎 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 N2. 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, N2, N2O, N2O4, NO, NO2, O, O2, and O3 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.

870631.fig.0015
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 reactionH2+12O2H2O,(68) while the IGE model solves the equilibrium reaction with nine product species,H2+12O2H,H2,H2O,H2O2,HO2,O,O2,O3,OH.(69) 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.

870631.fig.0016
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−1 K−1
Δ𝐻𝑓: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 =𝑚𝑗=1𝑛𝑗
(𝑔):Species in the gas phase
(𝑙):Species in the liquid phase
𝐻:Enthalpy, kJ
:Specific molar enthalpy of the reactant mixture, J mol−1
𝑗:Standard state specific molar enthalpy of species 𝑗, J mol−1
𝑇:Temperature, K
𝑆:Entropy, KJ·K−1
𝑠:Specific molar entropy of the reactant mixture, J mol−1 K−1
𝑠𝑗:Specific molar enthalpy of species 𝑗, J mol−1 K−1
𝑠𝑗:Standard state specific molar entropy of species 𝑗, J mol−1 K−1
𝑈:Internal energy, kJ
𝑃:Pressure, atm
𝑃𝑗:Partial pressure of species 𝑗, atm
𝑃:Standard state pressure = 1 atm
𝑅:Ideal gas constant = 8.314472 J·K−1·mol−1
𝑉:Volume,m3
𝑦𝑗:Mole fraction of species 𝑗
𝜂𝐶:Compressor isentropic efficiency
𝜙:Equivalence ratio
𝜇𝑗:Chemical potential of the 𝑗th species, J·mol−1
𝜇𝑗:Standard state chemical potential of the 𝑗th species, J·mol−1.

Acknowledgment

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

References

  1. P. Atkins and J. de Paula, Physical Chemistry, W.H. Freeman, New York, NY, USA, 8th edition, 2006.
  2. O. Redlich and J. N. S. Kwong, “On the thermodynamics of solutions,” Chemical Reviews, vol. 44, no. 1, pp. 233–244, 1949. View at Scopus
  3. 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. View at Scopus
  4. D. M. Kerrick and G. K. Jacobs, “A modified Redlich-Kwong equation for H2O, CO2, and H2O-CO2 mixtures at elevated pressures and temperatures,” American Journal of Science, vol. 281, no. 6, pp. 735–767, 1981. View at Scopus
  5. 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. View at Scopus
  6. 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.
  7. 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.
  8. 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.
  9. 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.
  10. F. J. Zeleznik and S. Gordon, “Calculation of complex chemical equilibria,” Industrial & Engineering Chemistry, vol. 60, no. 6, pp. 27–57, 1968.
  11. 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.
  12. 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. View at Scopus
  13. M. L. Michelsen, “The isothermal flash problem. Part I. Stability,” Fluid Phase Equilibria, vol. 9, no. 1, pp. 1–19, 1982. View at Scopus
  14. M. L. Michelsen, “The isothermal flash problem. Part II. Phase-split calculation,” Fluid Phase Equilibria, vol. 9, no. 1, pp. 21–40, 1982. View at Scopus
  15. 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. View at Publisher · View at Google Scholar · View at Scopus
  16. G. Soave, “Equilibrium constants from a modified Redlich-Kwong equation of state,” Chemical Engineering Science, vol. 27, no. 6, pp. 1197–1203, 1972. View at Scopus
  17. 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.
  18. 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.
  19. 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. View at Scopus
  20. 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. View at Scopus
  21. 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. View at Publisher · View at Google Scholar
  22. 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. View at Scopus
  23. S. I. Sandler, Chemical, Biochemical, and Engineering Thermodynamics, John Wiley & Sons, New York, NY, USA, 4th edition, 2006.
  24. 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. View at Scopus
  25. Y. A. Cengel and M. A. Boles, Thermodynamics—An Engineering Approach, McGraw–Hill, New York, NY, USA, 4th edition, 2002.
  26. 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. View at Scopus
  27. H. W. Cooper and J. C. Goldfrank, “B-W-R constants and new correlations,” Hydrocarbon Processing, vol. 46, no. 12, pp. 141–146, 1967.
  28. S. Bhattacharjee, The Expert System for Thermodynamics: A Visual Tour, Prentice Hall, New York, NY, USA, 2002.
  29. 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.
  30. 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. View at Scopus
  31. 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. View at Publisher · View at Google Scholar · View at Scopus
  32. 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. View at Scopus
  33. 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. View at Publisher · View at Google Scholar · View at Scopus
  34. 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.
  35. C. Paolini and S. Bhattacharjee, “Solving chemical equilibrium problems online,” Journal of Chemical Education, vol. 87, no. 4, p. 456, 2010. View at Publisher · View at Google Scholar
  36. 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.
  37. 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.
  38. 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.
  39. 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.
  40. V. Babrauskas, Temperatures in Flames and Fires, Fire Science and Technology, Washington, DC, USA, 2006.