Python port of MATLAB code used in MAE 250F
We are asked to start from the partition functions to compute equilibrium thermodynamic properties for air as a function of pressure and temperature. We will assume that the air is a mixture of N 2 , and O 2 at low temperature (using the approximate mole ratio of 79% to 21%). At high temperature, we can consider the air is a mixture of the following species: N 2 , O 2 , NO, N, O, N + , O + , and e − . The objective is to write a computer program to compute mixture’s specific enthalpy, h, internal energy, e, and entropy, s, as well as the mass fractions of the species, c i , compressibility factor, Z, and total density, ρ. In addition we shall find the specific heat ratio, γ, and speed of sound, a, of the gas mixture as functions of pressure, p, and temperature, T.
First the code will be validated through direct comparison of computed results with established results from literature. Following validation, more results will be given to prompt discussion of the thermodynamic properties at high temperatures.
For a gas in equilibrium, we can apply the Law of Mass Action for every reaction. K p (T) = ∏ i P ν i i = ( kT V ) ∑ i ν i ∏ i Q ν i i e −θ d /T (1) where P i is the partial pressure of species i, ν i is the stochiometric coefficient of species i, k is the Stefan-Boltzmann constant, V is the volume of the gas mixture, T is the equilibrium temperature, Q i is the partition function of species i, and θ d is the characteristic temperature associated with a dissociating gas. In the treatment of ionizing gases, we can simply replace θ d in Eq. (1) with θ I .
Equation (1) is the fundamental relationship that makes the entire study possible. The Law of Mass Action is so powerful for thermodynamics because it allow us to relate the partition functions, derived from a statistical treatment of the quantum energy levels, to the thermodynamic properties of partial pressure. The partition functions are solely functions of temperature and can therefore be calculated with ease. After finding the equilibrium constants, they can be related to the partial pressures and, along with some other equations, allow us to find all thermodynamic properties.
For our gas mixture, we will have five reactions that occur at different temperatures. They are N 2 . . . 2 N O 2 . . . 2 O NO. . . N + O O. . . O + + e – N. . . N + + e –
and the associated equilibrium constants (from Eq. (1)) are therefore
K p,N 2 = K 1 = P N 2 (P N ) 2 (2a) K p,O 2 = K 2 = P O 2 (P O ) 2 (2b) K p,NO = K 3 = P NO P N P O (2c) K p,O + = K 4 = P O P O + P e − (2d) K p,N + = K 5 = P N P N + P e − (2e)
From Dalton’s Law of partial pressures, we also know that the sum of all these partial pressures will equal the specified pressure of the system
P = P N 2 + P N + P O 2 + P O + P NO + P e − + P O + + P N + (3)
Combining Eq. (2) and Eq. (3) we have a system of 6 nonlinear equations for 8 unknown partial pressures. To find a complete set of equations we then simply need to apply conservation of charge which can be expressed in terms of pressure as
P e − = P O + + P N + (4) The second equation comes from a treatment of conservation of species which results in a form relating partial pressures as
N O N N = 2P O 2 + P O + P NO + P O + 2P N 2 + P N + P NO + P N + (5)
where we can find
N O N N = 0.21 0.79
for air at room temperature. With the addition of these two equations we now have a complete set of 8 nonlinear equations for 8 unknowns. The derivations of these equations are fairly straightforward but the details have been left out for the sake of brevity. The class notes of Prof. Zhong or the text book of Anderson are excellent resources for the details. Now it simply falls to requiring a calculation of the total partition functions for each species in the gas (shown next) which will provide the equilibrium constants needed on the left-hand-side of Eq. (2). Then we can use numeric tools to solve the system of equations. Once the partial pressures are calculated for each species at every temperature of interest, we can use gas relations to calculate all other thermodynamic properties.
The total partition functions come from multiplied-sums of the translational, rotational, vibrational, and electronic partition functions. From derivations based in statistics of quantum mechanics, these partition functions can be shown to be
Q t = V [ 2πmkT h 2 ] 3/2 (6a) Q r = T σθ r (6b) Q v = 1 1 − e −θ v T (6c) Q el = ∞ ∑ 0 g i e −θ el,i /T (6d)
Q = Q t Q r Q v Q el (7)
where the partition function for the electronic energy is found from n-term approximations for the infinite sum. For each species, the number of n-terms in the approximation varies by how influential the exponential term is at high temperatures. When monatomic molecules are considered, the rotational and vibrational partition functions are simply Q r = Q v = 1. For homo-nuclear diatomic molecules (N 2 and O 2 ), we have σ = 2 and for the hetero-nuclear diatomic molecule (NO), it is σ = 1.
The characteristic temperatures given in the partition functions, θ r , θ v , θ d , and θ I are given in Table 2. The degeneracy values (g i ) used are given in Table 3 and the associated electronic characteristic temperatures are listed in Table 4. Assuming I am capable of employing some numeric scheme to successfully calculate the partial pressure at each temperature, as mentioned we will use the partial pressures to find thermodynamic properties from those values. From pressure we can go to molar concentration as well as mass density directly
[x i ] = P i RT kmol/m 3 (8)
ρ i = P i R i T kg/m 3 (9)
where R = 8314 J/kmol-K is the universal gas constant and R i is the gas constant (in J/kg-K) for each species
R i = R M i (10)
where M i is the molecular weight of the species. The molecular weights of all the species considered in this study are given in Table 1. The total density is found from summation of the partial densities of each species, i,
ρ = ∑ i ρ i
which, once calculated, then also allows us to find the mass fraction
c i = ρ i ρ (12)
and the last quantity we calculate is the mole fraction, w i , which is defined as
w i = c i M M i (13)
M = ∑ i c i M i
is the average molecular weight of the mixture, and M i is the molecular weight of species i. Both have units of kg/mole.
After calculating all of the molar and mass fractions and concentrations, we proceed to calculating macroscopic thermodynamic properties such as the compressibility factor, Z, total specific internal energy, e, total specific enthalpy, h, and total specific entropy, s. In this study, the compressibility factor is found from
Z = P ρR 0 T (14)
where R 0 is the gas constant of the mixture at room temperature and is found from the known molar mass of air at room temperature, M 0 = 28.97 kg/kmol. The compressibility factor is an easily calculated macroscopic property and the dimensionless form was therefore used to validate the code (Figure 1) as will be discussed more completely later. To calculate the specific internal energy of every species, we can do it piece by piece. This is because the internal energy is calculated from the sum of different modes of energy: translational, e t , rotational, e r , vibrational, e v , electronic, e el , and the correction from the zero state, e 0 . For any species, i, the internal energy is
e i = (e t ) i + (e r ) i + (e v ) i + (e el ) i + (e 0 ) i (15)
The simplest energy to calculate is the translational energy; it is the exact same form for all species:
e t = 3 2 RT (16)
Next simplest is the rotational energy. The only thing to note is that monatomic molecules have no rotational energy, therefore it is calculated as
(e r ) N = (e r ) O = (e r ) e − = (e r ) O + = (e r ) N + = 0 (17a) (e r ) N 2 = R N 2 T (17b) (e r ) O 2 = R O 2 T (17c) (e r ) NO = R NO T (17d)
Next in line is the vibrational energy. Again, any of the monatomic molecules have no vibrational energy. The vibrational energy for diatomic molecules is found from the assumption in quantum mechanics that the diatomic molecules are a perfect spring-mass system. The vibrational energy is calculated as
(e v ) N = (e v ) O = (e v ) e − = (e v ) O + = (e v ) N + = 0 (18a) (e v ) N 2 = (θ v ) N 2 T e (θ v ) N 2 /T − 1 R N 2 T (18b)
(e v ) O 2 = (θ v ) O 2 T e (θ v ) O 2 /T − 1 R O 2 T (18c)
(e v ) NO = (θ v ) NO T e (θ v ) NO /T − 1 R NO T (18d)
The vibrational characteristic temperatures for the species are given in Table 2 In this simulation, the two-term approximation of the electronic mode for internal energy was used for every species.
(e el ) i = R i θ 1 ( g 1 g 0 ) i e −θ 1,i /T 1 + ( g 1 g 0 ) i e −θ 1,i /T (19)
where the degeneracy values and characteristic temperatures are given in Table 3 and Table 4. This approximation may not be sufficient at high temperatures and validation of the method shall hopefully show if this is the case. When calculating the internal energy, we have to be very careful how we set the zero state energy. e 0 and how we treat the repercussions of that choice on the form of e 0 . In the calculations we will use the common reference state of zero energy to be the completely combined diatomic state. Because of this choice, the calculation of internal energy is slightly more complicated than if we had chosen a completely atomic state at T = 0 K. The zero energy state corrections for the various species are
(e 0 ) N 2 = (e 0 ) O 2 = (e 0 ) e − = 0 (20a) (e 0 ) N = R N 2 (θ d ) N 2 (20b)
(e 0 ) O = R O 2 (θ d ) O 2 (20c)
(e 0 ) NO = R NO ( −(θ d ) NO + 1 2 [(θ d ) N 2 + (θ d ) O 2 ] ) (20d)
(e 0 ) N + = R N + 2 [(θ d ) N 2 + (θ I ) N ] (20e)
(e 0 ) O + = R O + 2 [(θ d ) O 2 + (θ I ) O ] (20f)
where the N 2 and O 2 subscripts on θ d indicate that in Table 2, as is standard in literature, the characteristic temperature for disassociation is attached to the diatomic molecule but we are using it in the monatomic species equation.
Once all the different modes of internal energy are calculated, they can be summed for each species, i, according to Eq. (15).
Following internal energy, there is a trivial calculation to find enthalpy. For any species, i, the enthalpy is
h i = e i + R i T (21)
The entropy of each species in the gas mixture is also found from the partition functions given in Eq. (6) and derived from the famous Boltzmann relation to have
S = Nk [ ln Q N + 1 + T ∂ lnQ ∂T ] (22)
If that calculation is carried out for all the modes of energy: translational, rotational, vibrational, and electronic the results are as follows
s t = 5 2 RlnT − Rlnp + R ln ( 2πM h 2 ) 3/2 k 5/2 + 5 2 (23a)
s r = R [ ln ( T σθ r ) + 1 ] (23b)
s v = R [ −ln ( 1 − e −θ v /T ) + θ v /T e θ v /T − 1 ] (23c)
s el = R [ lng 0 + ln ( 1 + g 1 g 0 e −θ 1 /T ) + ( g 1 g 0 ) ( θ 1 T ) e −θ 1 /T 1 + g 1 g 0 e −θ 1 /T ] (23d)
where the entropy is found for each species (the subscript i has been omitted above but should be imagined to be on each gas-specific term). In the calculation of the translational mode we have the mass per molecule M i . To differentiate this from the molecular weight, the script M, M i , will be used on a molar basis but when referring to individual particles/molecules/elements and the standard M, M i , is used. The mass of particle for each species is found the molecular weight and Avagadro’s number, N A = 6.022 × 10 23 particle/mol.
M i = M i N A (24)
Then the entropy for each gas species, i is found simply as
s i = (s t ) i + (s r ) i + (s v ) i + (s el ) i (25)
From the specific internal energies, enthalpies, and entropies of each species, the total energies can be found simply from applying mass-weighted sums
e = ∑ i c i e i (26a) h = ∑ i c i h i (26b) s = ∑ i c i s i (26c)
Thus we conclude our equations needed for the macroscopic thermodynamic properties required to be calculated for the study. The following will summarize the equations required to calculate the other various thermodynamic quantities studied.
In the course of the investigation, other thermodynamic properties were calculated for presentation; the specific heat at constant volume, c v , the specific heat ratio, γ, and the speed of sound, a. The specific heat was calculated from the partition functions (and a two-term approximation was used for the electronic mode contribution). For every species, the specific heat is built up from the quantum energy modes of translation, rotation, vibration, electronic
(c v ) i = ((c v ) t + (c v ) r + (c v ) v + (c v ) el ) i (27)
where for each species, i, we find
(c v ) t = 3 2 R (28a)
(c v ) r = R [ θ v 2T sinh θ v 2T ] 2 (28b)
(c v ) v = R (28c)
(c v ) el = R ( θ 1 T ) 2 ( g 1 g 0 ) e −θ 1 /T [ 1 + ( g 1 g 0 ) e −θ 1 /T ] 2
(28d) where the specific heat contribution to all the monatomic species from vibration and rotation are zero, much like the rotational and vibrational contributions to internal energy seen in Eq. (18) and Eq. (17), but the exact form for each species is not written out explicitly this time. Then the total specific heat is summed over the number of species, i
c v = ∑ i c i (c v ) i (29)
NOTE: This specific heat is for nonreacting gas mixtures. It was used incorrectly, see discussion on specific heat and specific heat ratio in the Discussion & Conclusion section. The specific heat ratio was found next as
c p = c v + ˆ R (30)
γ = c p c v (31)
where the mass-averaged gas constant is
ˆ R = ∑ i c i R i (32)
Once the specific heat ratio was calculated, finding the speed of sound was trivially
a = √ γ ˆ RT (33)
In reporting results, it is often more illustrative to present variables that have been normalized or non-dimensionalized in some way. Further, for the effort of validation, several non-dimensional values were calculated which will be described later, these are non-dimensional internal energy, e ∗ and non-dimensional entropy, s ∗ . The internal energy is
e ∗ = Ze ˆ RT (34)
and the entropy is
s ∗ = Zs ˆ R (35)
A non-dimensional enthalpy will also be reported which is similarly
h ∗ = Zh ˆ RT (36)
The species densities (from Eq. (9)) are normalized against the density of the the gas mixture (air) at room temperature (ρ 0 = P/[R 0 T 300 K ]).
ρ ∗ = ρ i ρ 0 (37)
The specific heat will also be nondimensionalized as
c ∗ v = Zc v ˆ R (38)
And lastly, the speed of sound is normalized against the ideal speed of sound a ideal =
√ γ 0 R 0 T (39) where γ 0 = 1.4. So a ∗ = a a ideal (40)
We have established here the equations that will be the framework for the code used to compute all desired thermodynamic properties. Next we can get into the algorithm and walk through the approach that was used in the program of this study.
The computer code was constructed to loop through both pressure and temperature as independent variables. The pressure was chosen at discrete values, chosen specifically so validation with Refs. [1, 2, 3] was possible. These were
P in = 0.0001 0.0010 0.0100 0.1000 1.0000 10.000 100.00
Temperature, on the other hand, was chosen on the domain T ∈ [2000,15000] with much smaller step sizes. The exact size of the step depended simply on how many steps were chosen, typically N = 1000.
The program ran by first looping through pressure. Then at each pressure, another for-loop would march through temperatures. Then at each temperature, all of the calculations discussed below are carried out. Once they are complete, the temperature loop takes another step and repeats the process. Of course, once the entire temperature span has been covered, the pressure is incremented to the next value and the entire temperature loop is carried out once more. In this way, the program automatically generates a heap of data that is used to validate and elucidate. Once the temperature is specified in the program, the partition functions for each mode of energy (given in Eq. (6)) are calculated for each species. The total partition function is then calculated for each species, as dictated by Eq. (7). With calculated partition functions, we find the equilibrium constants given in Eq. (1). As I was not specific initially on what stoichiometric coefficients, ν i , went into the Law of Mass Action, for clarity I will show the exact equations used in the program here
K 1 = 1 kT e [ (θ d ) N 2 /T ] Q N 2 (Q N ) 2 (41)
K 2 = 1 kT e [ (θ d ) O 2 /T ] Q O 2 (Q O ) 2 (42)
K 3 = 1 kT e [(θ d ) NO /T] Q NO Q N Q O (43)
K 4 = 1 kT e [(θ I ) O /T] Q O Q O + Q e − (44)
K 5 = 1 kT e [(θ I ) N /T] Q N Q N + Q e − (45)
and it is noted that these are consistent with the way the Law of Mass Action is written in terms of pressure in Eq.(2) which is vitally important to get correct results.
Now that we have the equilibrium constants, we simply must employ the 8 equations, described in the discussion above, for the 8 unknown partial pressures (Eqs. (2),(3),(4),(5)). To do this, we first create a matrix of partial pressures. The order of these is
p 1 → N 2 p 2 → N p 3 → O 2 p 4 → O p 5 → NO p 6 → e − p 7 → O + p 8 → N +
Then it can be simply shown that Eqs. (2),(3),(4),(5) are re-arranged as
0 = K 1 ( p 2 2 ) − p 1 0 = K 2 ( p 2 4 ) − p 2 0 = K 3 (p 2 p 4 ) − p 5 0 = K 4 (p 6 p 7 ) − p 4 0 = K 5 (p 6 p 8 ) − p 2 0 = 21 79 (2p 1 + p 2 + p 5 + p 8 ) − (2p 2 + p 4 + p 5 + p 7 ) 0 = p in − (p 1 + p 2 + p 3 + p 4 + p 5 + p 6 + p 7 + p 8 ) 0 = p 6 − (p 7 + p 8 )
When written in this form, the solution to the pressure matrix is easily solved with MATLAB’s built-in fsolve routine. The form of the system of equations written above was chosen to avoid any instances where a very large number was being subtracted from a very small number (or vice versa) that might introduce round-off or any other sort of numeric error into the solution. As it is, when the equilibrium constants K are very large, they will be balanced by being multiplied by a very small partial pressure.
To execute fsolve, an initial guess is necessary. For this study, it was always sufficient to provide a very first guess that the total pressure is shared equally by the total number of species (P i = P/I). fsolve uses this very first guess and generally had no problem providing a solution for the real partial pressures. Following that very first guess for the program, the solution given by fsolve was assigned to the guess of the next step. Because the partial pressure did not change dramatically between each step, the results of the last step are generally very similar to those for the next and fsolve could rapidly converge to solutions at each temperature step. This seemed to be an efficient method of solving the nonlinear equations over the span of temperatures.
Once the partial pressures are calculated, finding all other thermodynamic properties is straight- forward. The species densities and molar concentration are found immediately according to Eq. (8), Eq. (9). The total density is also found straightaway from a summation of those species densities. The mass fraction is calculated from the densities following Eq. (12). The mass fraction then allows us to calculate the molar fraction according to Eq. (13). A plot of the mole fraction with temperature is given in Ref.  and was thusly chosen as a plot in the results.
Following the calculation of the partial pressures and derivative values (mole concentration, density, etc.), the program runs through calculations of compressibility, internal energy, enthalpy, and entropy. The equations in the program follow directly those given for internal energy: Eq. (15), enthalpy: Eq. (21), and entropy: Eq. (23). The total values are found from Eq. (26).
The chosen programming interface for this study was MATLAB and the code is given in the Appendix. To compare the results with those from literature, WinDig was used to digitize plots off of PDF files. From WinDig, .DAT files were created that were able to be called by MATLAB in code and plotted to provide direct comparisons with the results of the simulation. Those results will be summarized next.
Before discussing the merits of the results output from the code, we must first validate that the solutions acquired are reliable and accurate. To that end, the results of the simulation are initially validated against those of Hansen from 1958. The first validation is done with the compressibility factor. This is a useful variable to compare as it can be computed simply from the thermodynamic properties of total density, the gas constant of air at room temperature, and the specified pressure and temperature. Three of those values are known so only the total density is dependent upon the solution of the system of equations defined above. That fact helps pinpoint errors, if there are any. The results of the validation are given in Figure 1. In this figure, we see excellent agreement at all temperatures and for all pressures. There is virtually no disagreement between the digitized results from Ref.  and the results of the current code. Any slight deviation can be assessed to the inaccuracy of digitalization process. This first result gives us confidence to continue validation.
After trusting that the code was producing reliable partial pressures and associated fluid properties, we proceed to calculate the specific, total internal energy, enthalpy, and entropy. Hansen reports a non-dimensional internal energy and non-dimensional entropy so for validation purposes, those values were calculated by the program in this study. The non-dimensional energy is given in Eq. (34). The first results of non-dimensional internal energy are given in Figure 2. The nondimensional energy tracks very similar paths at low temperatures but begin to deviate as the temperature increases. Because the compressibility was so strongly validated, I have no reason to think that the gas properties going into the internal energy calculation are incorrect. Instead, there must be some source of error in the calculation of internal energy itself.
I first suspected that this deviation is due to the two-term approximation to the electronic energy for internal energy. As a reminder, in calculating the partition functions (See Eq (6) – which eventually led us to the partial pressures at equilibrium temperatures) we used, in some cases, a 6 term approximation for the electronic mode. This is was done because the exponential terms at very high temperatures were not negligible for some species until that term. Therefore, such results as compressibility can be highly accurate while the internal energy, using a two-term approximation (Eq. (19)) is not. It is known that electronic energy has almost no contribution at low temperatures (this was shown in HW 4 of this class when similar electronic contribution to entropy at 1000 K was less than 4 %) but must become more influential to the total as temperatures rise. To test this, I compared the electronic mode of energy to the translational mode as temperatures and pressures varied. I created a variable e comp to compare these values: e comp = ∑ i c i (e el ) i (e t ) i (46)
where the translational mode alone was chosen because it is the only mode of energy that every species in the gas has and therefore served as a useful baseline for this test. This comparison variable was mass-weighted in the same way as the standard internal energy. The result of this calculation is plotted in Figure 3 and to me are very interesting. What we see in Figure 3 is a strong influence from pressure as to at what temperature the electronic energy peaks, dips down, and then rises again. Furthermore, that at higher pressures, the electronic energy rises to a value that is almost 25% of the translational energy mode. This indicates that the electronic energy mode becomes more important as temperature increases and implies that any error in the two-term approximation will be amplified at high temperatures. With this knowledge, we can review Figure 2 again and see that as pressure increases, the deviation of the internal energies increases (with increasing temperature); the error at p = 100 atm is much greater than p = 0.1 atm for T > 10000 which is exactly what Figure 3 would lead one to suspect. In that case, the code is not necessarily wrong, but simply requires a more precise method of calculating the electronic mode’s contribution to internal energy. Following the internal energy, entropy was validated against Ref. . The results are given in Figure 4. The same trends are seen here as in Figure 2; the error between the value calculated with the two-term approximation begins to deviate from the expected value at increased temperatures.
The validation efforts here have shown how the electronic mode of quantum energy can contribute to the proper calculation of properties. In the compressibility factor (Figure 1), where a 6-term approximation for the electronic mode was used, the results of the current study match extremely well with those from literature. Then in the internal energy and entropy calculations (Figure 2, Figure 4), where only a 2-term approximation was used, the results of the current study deviate from those of literature as temperature increases. Assuming that the lack of complete calculations of the electronic mode are the only source of error in the code, we can say that the algorithm at least is validated and can continue with the investigation of other fluid properties. In doing so, we expect that in any investigation based on the partial pressure calculations (from the 6-term approximation) to be highly accurate, while any reported values derived separately (with 2-term approximations) to be incorrect at high temperatures.
Figure 1: testing the image directory (delete this / move it to the correct location later) The first properties to report are the species compositions as the equilibrium temperature increases from 2000 K to 15000 K and as we vary the system pressure. The compositions will be reported in terms of mole fraction, w i (Eq. (13)), mass fraction, c i (Eq. (12)), and partial density, ρ i (Eq. (9)). These are reported in Figure 5, Figure 6, and Figure 7, respectively.
Discussion of the trends seen in the mole fraction, mass fraction, and density figures will bereserved for the following section. Until then, we turn our attention to other thermodynamic properties of interest. First, we view the enthalpy in Figure 8. We must note that because enthalpy was derived from internal energy, even without validation we can expect the results to be slightly inaccurate at high temperatures.
Following the calculation of enthalpy, we can see the results of the dimensionless specific heat in Figure 9. In the specific heat calculation, Eq. (28), we note that here too we used a two-term approximation for the electronic mode. As such we also expect the specific heat calculation to have some error at increased temperatures. The next interesting variable, based on the specific heat calculations given above, is the ratio of specific heats which is shown in Figure 10. Note again: the specific heat ratio was derived directly from c v so any error due to the two-term approximation of the electronic mode will also be apparent in γ.
Finally, we report the nondimensional speed of sound in Figure 12. The nondimensional value ranges from approximately unity at low temperatures up to nearly 2.2 at increased temperatures. For the last time: this value was also derived from two-term approximations on the electronic mode so it too will have some error at elevated temperatures.
As Hansen notes, because the compressibility is not influenced by vibrational excitation it has a value of unity until the moment oxygen begins to dissociate. Compressibility, in a sense, is also an indicator of a gas’s deviation from “ideal” and having a value of unity means the ideal gas law is applicable up to the temperature that dissociation begins. Oxygen is the first to dissociate and, having a molar ratio in air of approximately 21%, we expect and indeed see another plateau in Figure 1 at Z = 1.2. In our simplification of the air species we have only nitrogen and oxygen, and for this reason the compressibility increases to a value of 2.0 until nitrogen also completely dissociates. At this point, all diatomic molecules have dissociated into single atoms (for this discussion I am neglecting the production of any NO). As temperature further increases, the single atoms begin to ionize. Single ionization of the atoms doubles the number of gas particles under consideration again (one N produces an electron an N ion), so that eventually the compressibility approaches 4.0 when these reactions are complete.
The effect of pressure on the compressibility is to slow or accelerate the reaction processes indirectly through a shift in the dissociation properties. At higher pressures, the dissociation is delayed until higher temperatures (See, for instance Figure 5). In Figure 1 we see for a system pressure of p = 0.0001 atm, the dissociation and ionization is nearly complete at T = 11000 K. However, at a system pressure of p = 100 atm, at T = 14000 K we have only just seen a dissociation of N 2 into N.
The other interesting feature to discuss from the compressibility chart is spurred from discussion in . The effects of chemical reactions are most intense where the slope of the compressibility is a maximum, or at the point of inflection in the compressibility slopes. We can compare the temperature at points of inflection of Figure 1 for compressibility to match with temperatures in Figure 5 at which we see the sharpest changes in composition. Also interesting to note that the slope of Z is nearly zero at the transition from one major reaction to another. This shows that one reaction is essentially complete before the next reaction starts. For the highest pressures, the slope of Z never quite flattens to zero and indeed we see in Figure 5 that for p = 100 atm that N 2 begins to dissociate to N before the O 2 to O reaction is complete but even then the overlap is small. This observation solidifies the assumptions made in some analysis of gas mixtures that the reactions are all occurring independent of the other species in the mixture. The internal energy, as given in Figure 2, shows the same trend in dissociating gas species as that seen in the familiar change of internal energy of H 2 O as it changes phases. Included for reference is a chart from Ref.  (Figure 13). The difference between a phase change and dissociation is the constant temperature of phase change that does not exist for dissociation. However, the trends are the same. When nitrogen dissociates from a diatomic molecule into a monatomic element for instance (between 4000 K and 5000 K for p = 0.0001 atm on Figure 2) we have only a 25% increase in temperature while experiencing a 300% increase in the internal energy.
It is this very energy absorption property of dissociating (and ionizing) gas species that differentiates the results of thermally perfect calculations from calorically perfect calculations. Like Prof. Zhong mentioned in class, calorically perfect models of gas predict very high temperatures across shocks when in reality the temperatures are much lower. The lower temperature is thanks, in part, to the energy absorption properties of dissociation.
The increase in entropy between steps of dissociation seen in Figure 4 can be thought of as the increase in disorder associated with a doubling of the number of particles in the system. Increasing temperature results in mild increases the entropy as the accountability of any single particle in the mixture becomes harder to track. But for every subsequent dissociation of one particle into two, or ionization of one particle into two, we have very large increases in the reported entropy value. The results match our expectations of entropy; the value never decreases for increased temperature.
The results of species concentrations (molar, mass, mass per volume) all indicate the relative contributions of species to the net thermodynamic properties. The results shown in Figure 5, Figure 6, and Figure 7 are discussed, essentially, during the discussion of compressibility. It is reiterated here that the effect of pressure is to delay the dissociation and ionization reactions. For very low pressures (p = 0.0001 atm) complete ionization occurs just after T = 10000 K. By comparison, for p = 100 atm, the ionization of nitrogen is barely beginning by T = 15000 K. This result is expressed in all three of the relative species composition figures. Next, focusing our attention on p=1 atm, we see that diatomic oxygen begins dissociation at approximately 2000 K and is nearly complete by 4000 K. The presence of NO is rather insignificant as it only appears in trace amounts and flits in and out of existence between 3000 and 5000 K. Nitrogen dissociation completes around 9000 K at which point there is a small window of temperature of chemical equilibrium before ionization of oxygen and nitrogen begin to advance. This indicates that ionization can be ignored if a system at 1 atm pressure will experience under 10000 K temperature environments and that dissociation can be ignored for a system under 2000 K temperature environments.
The results of enthalpy follow from the discussion of internal energy and all major trends of internal energy are also present in Figure 8. The only difference between internal energy and enthalpy is the addition of the pressure term indicating how much work the system must perform to maintain the given volume through displacement of the environment.
The last plot of interest is the specific heat ratio, γ, shown in Figure 10. For air at room temperature it is a given that γ 0 = 1.4. Using the calorically perfect model of gas, this is a constant value for all temperatures. This number comes from the fact that air is mostly diatomic nitrogen and that the specific heat of a diatomic molecule that is fully excited translationally and rotationally is
c v = 5 2 , which means c p = 7 2
and thus γ = 1.4. As temperature increases for a diatomic molecule, the vibrational mode becomes fully excited and as the limit of T → ∞ we have
c v = 7 2 thus c p = 9 2
and γ = 1.29 – ignoring electronic mode-contributions and assuming the gas does not dissociate or ionize. However, for a monatomic gas, there is no rotational mode so if the translational state is fully excited,
c v = 3 2 and c p = 5 2
so γ = 1.667 for all temperatures (there is no vibrational mode either). From these two idealized situations – solely diatomic molecules and solely monatomic molecules – we have limits on what we expect for γ. It should be about γ = 1.29 on the low end (when the species is dominated by diatomic molecules), and approach γ = 1.667 when complete dissociation to monatomic molecules occurs. A tabular summary of expected values is presented in Table 5.
Amazingly (or perhaps not so amazing), these are nearly the precise limits seen in Figure 10. At T = 2000 K the vibrational mode of the diatomic molecules is very excited and no dissociation has yet occurred. Therefore, γ should be approaching 1.29 from above and, indeed, for all pressures we see γ ≈ 1.3. As temperature increases and dissociation splits all the nitrogen and oxygen into monatomic elements, the weighting of γ quickly increases towards the 1.667 limit. Let us follow the green line of p=1 atm. We see a spike (γ ≈ 1.47) around T = 7000 K which is right in the middle of the steepest dissocation rate of nitrogen. As temperature increases further, the specific heat continues on the trend toward a maximum value of γ ≈ 1.63.
It was pointed out to me just after writing this report that the equations used to build the specific heat are not accurate for a changing gas mixture. For one thing, the total specific heat written in Eq. (29) can be re-written as
c v = ∑ i c i (c v ) i = ∑ i [ c i ( ∂e i ∂T ) V ] (47)
however, the total specific heat can also be written as c v = ( ∂e ∂T ) V = ∂ ∂T ( ∑ i c i e i ) V = ∑ i [ c i ( ∂e i ∂T ) V + e i ( ∂c i ∂T ) V ] (48)
so we can compare the right hand sides of the equations in Eq. (47) that was used in the program with that of Eq. (48). Clearly there is a difference of a factor
e i ( ∂c i ∂T ) V
inside the summation. This term takes into account that the mass fraction changes with temperature along with the internal energy. Being a derivative of mass fraction, that term is completely uninfluential to the total specific heat at times when every species mass fraction is constant in temperature (at points of plateaus in Figure 6) when the change of mass fraction is zero and likewise most influential during large changes in c i seen in Figure 6. Therefore, for p = 1 atm at T = 7000 K (as pointed out above) the specific heat ratio will be most erroneous. Thus no physical interpretation of these regions should be made, but at the points of plataeu in Figure 6, the specific heat ratio should still be mostly accurate. The difference between Eq. (29) and Eq. (48) is not that the former is incorrect, it’s that it was applied incorrectly. Equation (29) is only for nonreacting gas mixtures. Whereas for reacting gas mixtures I should have used Eq. (48).
I also did some numeric experimenting to check the influence of the electronic mode on the specific heat ratio (which was neglected in the discussion of ideal limits shown in Table 5) and created a comparison. For the comparison, I created a plot of the specific heat ratio when no electronic mode was included. This comparison is given in Figure 11. In this plot we see the dashed lines of γ go straight up to γ = 1.667 which is the exact limit predicted for solely monatomic elements (which is essentially all we have at the high temperatures). This proves the expected theory. As for the contribution of the electric mode, I have no qualitative discussion other than to say it must be to blame for the behavior between what is shown in dashed and solid lines in the figure.
Also included is the normalized speed of sound. The dimensionless value is approximately unity at T = 2000 K and has multiple plateaus as first oxygen completely dissociates, then nitrogen does the same, and finally when the ionization of the atomic species is complete. Ultimately, the speed of sound for the completely ionized species is nearly 2.2 times greater than the speed of sound at standard conditions (and T = 300 K).
This study went into great detail explaining how thermodynamic properties can be attained for high temperature gases by approaching them from a statistical thermodynamics point of view. With partition functions calculated (as only functions of temperature), partial pressures of species in a gas mixture are (relatively) easily calculated. From partial pressures at equilibrium temperatures, all sorts of derivative thermodynamic calculations are possible. The numeric algorithm developed in this study was highly accurate at calculating partial pressures. However, the two-term electronic mode approximations on energy, entropy, specific heat - and the associated derived terms - suffered inaccuracies at high temperatures when the electronic mode was not negligible compared to the other modes. All the same, many interesting observations were capable from the results of this study.
One last note concerns the species reactions involved in the study. The program had been constructed with an a-priori knowledge of the mixture reactions. This knowledge allowed us to construct the species equations and associated equilibrium constants (See, for example Eq. (2)). For air, this is acceptable because of how widely studied the mixture has been and it is understood what reactions are producing what species in the gas. However, for some other unknown composition, it may be necessary to include a larger number of equilibrium constants (based on other possible reactions) and then allow the math to indicate whether or not the reactions are important.