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

Authors: Christopher P. Paolini; Subrata Bhattacharjee
Journal: Journal of Thermodynamics
 (2012)

---

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


---

## Body

## 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)pv=RT.
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)pv=zRT,
with z 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 z when computing the state of a particular gas.The ideal gas law was first introduced by Clapeyron in his 1834 work entitledMé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 v the volume b where b 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 a representing the constant of proportionality, van der Waals’ model is given by(3)P=RTv-b-av2.
For a given species, the constants a and b 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(4)(∂P∂v)Tc=0,(∂2P∂v2)Tc=0,
one will obtain two equations in the two unknowns, a and b, which can then be computed from knowledge of the species’ critical temperature and pressure, Tc,  Pc, 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)P=RTv-b-aTv(v+b).
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 a and b 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 Z 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)P=RTv-b-aαv(v+b),
includes an additional factor α that is multiplied by the attractive parameter a 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)P=RTv-b-aαv(v+b)+b(v-b),
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 Tr with κ being a function of the acentric factor ω. Optimal values for κ1, κ2, and κ3 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 flashmethod 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(8)P=RuTv̅2(1-cv̅T3)(v̅+b)-Av̅2,A=A0(1-av̅),B=B0(1-bv̅),
where A and B are empirical parameters, it is fairly accurate for species with densities up to 0.8ρcr [25]. The Benedict-Webb-Rubin EOS [26],(9)P=RuTv̅+(B0RuT-A0-C0T2)1v̅2+bRuT-av̅3+aαv̅6+cv̅3T2(1-γv̅2)e-γ/v̅2,
where the constants a, A0, b, B0, c, C0, α, 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 [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 pv=RT 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., cp=cv=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 p1=1atm and T1=3000K. The notation p1 and T1 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.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 Figure3 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 O2molecule 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 O2→O+O is defined as the standard state enthalpy change for the reaction at a specified temperature,(10)D∘=2ΔHf∘(O)-ΔHf∘(O2).
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.(11)O2⇌2O.
The amount of O2 molecules and O atoms present in equilibrium at a particular temperature is defined by the partial pressure equilibrium constant Kp 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 Kp for reaction (11) in terms of the partial pressures of the reactant O2 and the product O,(12)Kp=pO2pO2=yO2(P∘)2yO2P∘=[2α/((1-α)+(2α))]2(1-α)/((1-α)+(2α))P∘=4α21-α2P∘.
Figure 4 shows a plot of the log of Kp 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 Kp is constant and independent of temperature. The negative slope indicates reaction (11) is endothermic. As the temperature increases, Kp 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 Kp.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.Figure 4
Log of the partial pressure equilibrium constantKp versus the reciprocal of temperature for reaction (11).Using the van’t Hoff equation,(13)[∂ln⁡Kp∂(1/T)]p=-ΔH∘R,
and the slope of the log Kp v. 1/T plot, we can solve for ΔH°, the standard state molar enthalpy of reaction, and investigate how ΔH° varies with temperature. Figure 5 shows a plot of ΔH° 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.Figure 5
Standard state molar enthalpy of reactionΔH° versus temperature for reaction (11). ΔH°>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 functionG is defined as(14)G=H-TS=(U+PV)-TS=∑j=1mμjnj.
A system will be in equilibrium when the change in Gibbs free energy of the system vanishes to zero.Differentiating (14) we obtain(15)dG=dU+PdV+VdP-TdS-SdT.
At a constant temperature T and pressure P, (15) reduces to(16)dG=dU+PdV-TdS.
The combined first and second laws of thermodynamics requires(17)dU-TdS+PdV≤0⟹dG≤0,
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 j,(18)μj=μj∘+RTln⁡(PjP∘).
Substituting (18) into (14) and dividing by RT to obtain a dimensionless form, we obtain(19)GRT=∑j=1m[njμj∘RT+njln⁡(PjP∘)],
but(20)ln⁡(PjP∘)=ln⁡(PjPPP∘)=ln⁡(yjPPPP∘)=ln⁡yj+ln⁡(PP∘),
where the last term, ln(P/P∘), reduces to a constant since the system total pressure P is given as an input to the problem. Using (20), (19) now becomes(21)GRT=∑j=1m[njμj∘RT+njln⁡yj+njln⁡(PP∘)].
The minimum stationary point of (14) will be the vector of species molar values n⃗where dGvanishes. Differentiating G in (14) we obtain(22)dG=∑j=1mnjdμj+∑j=1mμjdnj,
but from the isothermal, isobaric Gibbs-Duhem equation we know that(23)∑j=1mnjdμj=0,
and so we seek the unique vector n⃗such that(24)∑j=1m[μj∘RT+ln⁡nj-ln⁡n+ln⁡(PP∘)]dnj=0.
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(25)∑j=1mAi,jnj=bi,i=1,…,a.
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(26)nj≥0,∀j,1≤j≤m.
We can express the atomic constraint equations given by (25) and (26) as equality constraints(27)ϕi=ϕi(n1,n2,…,nm)=∑j=1mAi,jnj-bi=0.
There will be a equations of the form (27) for each atom i present in the system. We define the Lagrangian function  ℒ=ℒ(n⃗,λ⃗) as
(28)L(n⃗,λ⃗)=μ(n⃗)+∑i=1aλiϕi.
The solution of (24) is thus the minimum stationary point where the gradient of (28) vanishes, that is,(29)∇L(n⃗,λ⃗)=∑j=1m(∂μ∂nj)p,T,ni≠jdnj+∑i=1aλi∑j=1m∂ϕi∂njdnj=0⃗.
The TEST IGE model uses a Newton-Raphson technique to solve (29). By expanding (29) we can write(30)∇L(n⃗,λ⃗)=∑j=1m(∂μ∂nj)p,T,ni≠jdnj+∑i=1aλi∑j=1m∂ϕi∂njdnj=∑j=1mμjdnj+∑i=1aλi[∂ϕi∂n1dn1+∂ϕi∂n2dn2+⋯+∂ϕi∂nmdnm]=∑j=1mμjdnj+λ1[∂ϕ1∂n1dn1+∂ϕ1∂n2dn2+⋯+∂ϕ1∂nmdnm]+λ2[∂ϕ2∂n1dn1+∂ϕ2∂n2dn2+⋯+∂ϕ2∂nmdnm]+⋯+λa[∂ϕa∂n1dn1+∂ϕa∂n2dn2+⋯+∂ϕa∂nmdnm]=[μ1+λ1∂ϕ1∂n1+λ2∂ϕ2∂n1+⋯+λa∂ϕa∂n1]dn1+[μ2+λ1∂ϕ1∂n2+λ2∂ϕ2∂n2+⋯+λa∂ϕa∂n2]dn2+⋯+[μm+λ1∂ϕ1∂nm+λ2∂ϕ2∂nm+⋯+λa∂ϕa∂nm]dnm=0⃗.
In order to satisfy (30), all the terms in the square brackets must vanish to 0. We therefore seek values of Lagrange multipliersλi⃗ such that(31)μj+∑i=1aλi∂ϕi∂nj=0,
for all j. By observing(32)λi∂ϕi∂nj=λiAi,j,
solving (30) amounts to solving a system of m equations for each species j where the jth species equation is given by(33)μj+∑i=1aλiAi,j=0,
and additionally, a population constraint equations for each element i where the ith constraint equation is given by(34)∑j=1mAi,jnj-bi=0,
and finally one constraint equation governing the total number of moles present in the system,(35)∑j=1mnj-n=0.
From (33), (34), and (35), we must solve a system of m+a+1equations in a+m+1unknowns. The IGE model uses the Newton-Raphson method to iteratively solve this system of equations. To construct our system matrix we let fj=fj(x⃗)=fi(x1,x2,…,xn) represent the jth species equation of (33), or(36)fj=μj+∑i=1aλiAi,j=0.
Expanding μj using (21) we have(37)fj=μj∘RT+ln⁡nj-ln⁡n+ln⁡(PP∘)+∑i=1aλiRTAi,j=0.
By letting(38)x1=ln⁡nj,x2=ln⁡n,x3=-λ1RT,x4=-λ2RT,…,xa+2=-λaRT,
we will obtain(39)∂fj∂x1=∂fj∂(ln⁡nj)=1,∂fj∂x2=∂fj∂(ln⁡n)=-1,∂fj∂x3=∂fj∂(-λ1/RT)=-A1,j,…,
and our Newton-Raphson equation becomes(40)∑i=1n∂f∂xiδxi=-f(x⃗)⟹δln⁡nj-δln⁡n+∑i=1aAi,jλiRT=-μjRT-∑i=1aλiRTAi,j.
We will have m equations of type (40) for each species j. For the population constraint equations given by (34), we let fi=fi(x⃗)=fi(x1,x2,…,xn) represent the ith population constraint equation such that(41)fi=∑j=1mAi,jnj-bi=∑j=1mAi,jexp⁡(ln⁡nj)-bi=0.
Then(42)∂fi∂xj=∂fi∂(ln⁡nj)=Ai,jnj,
and our Newton equation for (34) becomes(43)∑j=1mAi,jnjδln⁡nj=bi-∑j=1mAi,jnj.
We will have a 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 f=f(x⃗)=f(x1,x2,…,xn) represent (35),(44)f=∑j=1mexp⁡(ln⁡nj)-exp⁡(ln⁡n)=0,
where(45)∂f∂xj=∂f∂(ln⁡nj)=nj,∂f∂(ln⁡n)=-n,
and thus the Newton equation for the total system moles constraint becomes(46)∑i=1n∂f∂xiδxi=-f(x⃗)⟹∑j=1mnjδln⁡nj-nδln⁡n=exp⁡(ln⁡n)-∑j=1mexp⁡(ln⁡nj)=n-∑j=1mnj.
We can represent our system of Newton equations in matrix form by combining (40), (43), and (46) to obtain(47)[10⋯0-A1,1-A2,1⋯-Aa,1-101⋯0-A1,2-A2,2⋯-Aa,2-1⋮⋮00⋯1-A1,m-A2,m⋯-Aa,m-1A1,1n1A1,2n2⋯A1,mnm00⋯00A2,1n1A2,2n2⋯A2,mnm00⋯00⋮⋮Aa,1n1Aa,2n2⋯Aa,mnm00⋯00n1n2⋯nm00⋯0-n][δln⁡n1δln⁡n2⋮δln⁡nm-λ1RT-λ2RT⋮-λaRTδln⁡n]=[-μ1RT-∑i=1aλiAi,1RT-μ2RT-∑i=1aλiAi,2RT⋮-μmRT-∑i=1aλiAi,mRTb1-∑j=1mA1,jnjb2-∑j=1mA2,jnj⋮ba-∑j=1mAa,jnjn-∑j=1mnj].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 (“hp”) or constant pressure, constant entropy (“sp”) 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 (“tp”) problem until the specific enthalpy of the product mixture equals the specific enthalpy of the reactants. The direct hp constraint equation is derived as follows. Let h represent the dynamically changing specific enthalpy of the product mixture and h° represent the specific enthalpy of the reactants, a constant. Upon convergence of (47) we desire(48)h=h∘.
The specific enthalpy of the product mixture h is determined after each iteration of (47) through the summation(49)h=∑j=1mnjhj∘,
where hj∘is the standard state specific molar enthalpy of species j. As the number of moles of each species nj varies, the current value of the product mixture enthalpy varies as(50)∑j=1mnjhj∘RTΔln⁡nj.
For an hp problem, temperature Tvaries and approaches the adiabatic flame temperature as (47) converges. As T varies, the product mixture enthalpy also varies as(51)[∑j=1mnjcp,j∘R]Δln⁡T,
where cp,j∘ is the standard state specific molar heat capacity of species j. Equation (51) is clear sinceΔh=cpΔTfor an ideal gas. The standard state specific molar enthalpy and heat capacity in (50) and (51) are cast into dimensionless form by dividing by RT and R, respectively. Our additional constraint equation is thus(52)∑j=1mnjhj∘RTΔln⁡nj︸Enthalpychangethroughachangeinspeciesmoles+[∑j=1mnjcp,j∘R]Δln⁡T︸Enthalpychangethroughachangeintemperature+hRT︸Enthalpyoftheproductmixtureatthecurrentiteration=h∘RT︸Enthalpyofthereactants,aconstant
or(53)∑j=1mnjhj∘RTΔln⁡nj+[∑j=1mnjcp,j∘R]Δln⁡T=h∘-hRT.
Equation (53) is inserted into (47) and repeatedly solved yielding correction terms Δlnnj  andΔlnT. The final temperature T 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 sp problems. In this case we require the entropy of the product mixture to equal the entropy of the reactants, or(54)s=s∘,
where s represents the dynamically changing specific entropy of the product mixture and s° represents the specific entropy of the reactants, a constant. The specific entropy of the product mixture sis determined after each iteration of (47) through the summation(55)s=∑j=1mnjsj,
where sj is the specific molar enthalpy of speciesj given by(56)sj=sj∘-Rln⁡njn-Rln⁡P,
and sj∘ is the standard state specific molar entropy of species j given by(57)s∘=∫0Tcp(T)TdT.
As the number of moles of each species nj varies, the current value of the product mixture entropy varies as(58)∑j=1mnjsjRΔln⁡nj.
For an sp problem, temperature T varies and approaches an isentropic temperature as (47) converges. As T varies, the product mixture entropy also varies as(59)[∑j=1mnjcp,j∘R]Δln⁡T=∑j=1mnjΔln⁡Pj=∑j=1mnjΔln⁡njn+∑j=1mnjΔln⁡PP∘︸=0=∑j=1mnjΔln⁡nj-(∑j=1mnj)Δln⁡n=n-∑j=1mnj.
The last equality in (59) is from (46). Our additional constraint equation for an isentropic, isobaric process is thus(60)∑j=1mnjsjRΔln⁡nj+[∑j=1mnjcp,j∘R]Δln⁡T=s∘-sR+n-∑j=1mnj.
The second method used by the TEST IGE model to solve a hp or sp problem is an iterative method in which repeated constant pressure, constant temperature (i.e., tp) 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.Algorithm 1
Setα=0.4Specify the specific enthalpyh°  and pressure p of the reactantsCompute the temperature as a function ofh° and p,T=f(h°,p)Compute the equilibrium distribution of the product mixtureaccording toT and pRepeatSetT’=TCompute the thermodynamic state of this product mixturebased onh° and p giving a new temperature TSetT=αT+(1-α)T’Compute the equilibrium distribution of the productmixture according toT and pUntil|T’-T|<τwhereτ is a convergence error tolerance. NASA CEA [11, 36] uses the direct approach while STANJAN [39] performs repeated tp 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. Figure6 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, Figure6 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 labeledT2.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(61)CH4+2(O2+3.76N2)⟶CO2+2H2O+7.52N2.
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 setp1=100kPa and T1=298.15K. 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 h1.Once the specific enthalpy of the reactants is known, the IG model can be used to find the adiabatic flame temperature by definingState 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 p1=p2and h2=h1 to constrain pressure and specific enthalpy between the two states. Click the Calculate button again and the computed temperature displayed in the T2 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 theComposition 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 T2 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 theComposition Panel. The IGE model computes a lower adiabatic flame temperature than the IG model.All of the screen images of daemons shown in Figures1 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 Figure11. 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 constantKp of a reaction at different temperatures. IGE model-based daemons will report ln Kpin 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, ΔG°. From this value, ln Kp is easily found at a particular temperature T by(62)ln⁡Kp(T)=-ΔrxnG∘RT.
Figure 13 shows a plot comparing ln Kp versus the reciprocal of temperature for the equilibrium reaction(63)CO2⇌CO+12O2.
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,(64)ϕ=mfuel/moxidizer(mfuel/moxidizer)stoichiometric,
where m 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 is(65)C3H8+5(O2+3.76N2)⟶3CO2+4H2O+18.80N2.
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.Figure 12
TEST IGE model daemons will compute and report the log of the partial pressure equilibrium constantKp of a reaction in the upper right corner of the Composition Panel.Figure 13
Log of the partial pressure equilibrium constantKp   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(66)CH4+2O2⟶products,
where productsvaries 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 SpeciesIterative Method Wall Time (s)Direct Method Wall Time (s)Number oftp IterationsIterative Method Flame Temperature (K)Direct Method Flame Temperature (K)Difference2281.803.19305166516603322.4215.19305166516604216.173.3321384238420592.783.7083376337606455.584.474132623262072429.08†5.052973122312208*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 Table1.
Number of Product SpeciesSpecies Set2CO2, H2O3CO2, H2O, CO4CO2, H2O, CO, O25CO2, H2O, CO, O2, OH6CO2, H2O, CO, O2, OH, O7CO2, H2O, CO, O2, OH, O, H28CO2, H2O, CO, O2, OH, O, H2, H9CO2, H2O, CO, O2, OH, O, H2, H, O3Using the IGE model to perform asp 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 T2. For values of T2>Ts, where Ts is the exit temperature reached if the inlet air was compressed reversibly and adiabatically (i.e., isentropically), we can compute the compressor isentropic efficiency, ηC, using the IG and IGE model within the TEST Open System, Steady-State, Single Flow Equilibrium daemon.The compressor isentropic efficiency is defined as(67)ηC=IsentropiccompressorworkActualcompressorwork=wswa=ΔhsΔha=h2s-h1h2a-h1,
where h2s is the specific enthalpy of the exit state for an isentropic compression and h2a is the specific enthalpy of the exit state for an actual compression. The term Δhs 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 Δha 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 O2and 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 ηC 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, Figure16 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(68)H2+12O2⟶H2O,
while the IGE model solves the equilibrium reaction with nine product species,(69)H2+12O2⇌{H,H2,H2O,H2O2,HO2,O,O2,O3,OH}.
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 fromhttp://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 solvinghp and sp 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 tp 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].



---

*Source: 870631-2012-03-05.xml*