Thermal and Transport Properties for the Simulation of Direct-Fired sCO₂ Combustor

The fundamental research and technology development for direct-fired supercritical CO₂ (sCO₂) power plants is gaining the attention of researchers in academic institutions and industries, due to its remarkable theoretical promise of efficiency, compactness, and eco-friendliness. The current state-of-the-art optimized sCO₂ operating conditions (typical mixtures) are 200–300 atm pressure, 800–1000 K inlet temperature, and 1200–1500 K outlet temperature [1]. These unconventional and challenging operating conditions entail the need to review the fundamentals of all engineering fields involved in gas turbine power plant combustor technology. At these supercritical operating pressures, the mean free path between the molecules reaches to a distance where the intermolecular forces become prominent [2]. Though there always exists an attraction or repulsion force between the molecules of any fluid, at subcritical pressures the magnitudes of these forces are trivial [3] and do not impact the thermal properties of the fluid. In a supercritical combustion environment where there are hundreds of species and radicals, the influence of intermolecular forces depends on the proportion of the mixture constituents. Therefore, the rules or models of thermal properties which are developed to consider the combined effect of the mixture constituents are very important and need to be identified before a simulation. The thermodynamic behavior of supercritical fluids has been well recorded by researchers for applications such as petroleum, food processing, pharmaceutical, textile, metallurgical, and rockets. There is significant literature available on thermal properties, theoretical modeling, and advanced numerical simulations [4–10] of rocket combustion applications. However, for emerging sCO₂ combustion, the fundamental thermal quantities such as specific heats, the speed of sound, enthalpy, entropy, compression factor, etc., are not yet accurately quantified. Therefore, the main motivation of this work is to provide a quantification of the important thermal properties that are useful in calculations of combustion systems. The equation of state (EOS) is most important in quantifying the thermal properties of a system because every thermal property is directly related to pressure, temperature, and specific volume of that system, and the EOS is the relation that connects these three parameters. The first and foremost aspect of any supercritical combustion simulation is choosing a suitable EOS. Therefore, a brief review of EOSs and the formulations of EOS for combustion mixtures are discussed in detail, and then, a comprehensive analysis is carried out between most popular EOS such as Peng–Robinson (PRS) and Soave–Redlich–Kwong (SRK) at various turbulence levels to identify the best suitable EOS for sCO₂ combustion simulations.

In the Compression Factor (Z) section, the importance of compression factor (Z) in sCO₂ combustion modeling is discussed in detail, and a correlation to model Z in one-dimensional (1D) combustor analysis is proposed. Later, the important combustion parameters such as constant pressure specific heat (C_p), constant volume specific heat (C_v), ratio of specific heats $(γ)$ compression factor $(Z)$⁠, isothermal compressibility $(βT)$⁠, and process index $(ns)$ are quantified.

Finally, the high-pressure viscosity and thermal conductivity models for mixtures and pure components are reviewed from the literature, and the suitable models are identified for sCO₂ combustion simulations based on their accuracy and computational time.

It is known that, under supercritical conditions, the ideal gas assumption will not predict the system density correctly due to the existing repulsive or attractive forces between the molecules. These intermolecular forces completely alter the thermal and kinetic properties of combustion such as enthalpy, entropy, reaction rates, and so on; therefore, the kinetic tool which is being used for the simulation of sCO₂ combustion must have the capabilities to calculate the thermal and kinetic properties based on the EOS. The CHEMKIN-RG (extended version of CHEMKIN) is one such tool developed by Schmitt et al. [11] in 1993, and this tool is coupled with an in-house premixed conditional moment closure code (PCMC) [12] for estimating the state-related thermal properties of sCO2 combustion at different turbulence levels. The PCMC model as shown in Eq. (1), is a premixed turbulent combustion model which conditions the species mass fractions on the reaction progress variable (RPV) and closes the chemical source terms in the enthalpy equation with conditioned reaction rates [13]. In other words, the PCMC can estimate all the species mass fractions involved in combustion as the reaction progresses from unburnt to fully burnt conditions, at various turbulence levels by solving the second-order PCMC ordinary differential equation for each species with a two-point boundary value problem solver.

A high-pressure methane kinetic mechanism, Aramco 2.0 [14], which is validated up to 260 atm pressure in 67% argon dilution and up to 180 atm in 55% N₂ dilution, is used. The mechanism is also validated with 90% CO₂ diluted, 110 atm, unpublished shock tube auto-ignition data. The Aramco mechanism is also compared with the GRI 3.0 [15] and USC [16] mechanisms up to 40 bar with equivalence ratios up to 3 and found to be better performing against the available shock tube data [17]. This Aramco 2.0 mechanism, tailored for C1–C2 compounds has 73 species and 426 reactions.

Here, $ρ¯$ is the density, $Ñ$|ζ is the conditioned scalar dissipation (level of small-scale turbulence), $S̃c|ζ$ is the conditioned source term for the RPV equation, $ω̃˙i|ζ$ is the conditioned reaction rate, $Q̃i$ are the conditioned mass fractions, and the derivatives are with respect to the RPV [2,12].

The EOS is a correlation that describes the relation between pressure (P), temperature (T), and density (ρ) of a thermal system. Every fluid or mixture exists under a phase, i.e., solid, liquid, or gaseous vapor, based on the pressure and temperature acting on that system. At lower pressures and temperatures, the liquid phase of the system will have a different density than the gaseous phase. But, as the pressure and temperature increases, the difference between the phase densities decreases, and at one combination of pressure and temperature, the density of both phases becomes the same and at this point, the phase of the system cannot be distinguished based on the density. This combination of pressure and temperature is called the critical point, and any pressure and temperature combination above the critical point is considered as a supercritical state. A schematic diagram is shown in Fig. 1, to illustrate the critical point of CO₂ and the operating conditions of sCO₂ combustor in that state diagram. In these supercritical states, the molecules are pushed very close in such a way that the attractive or repulsive forces between the molecules become significant. A supercritical fluid has some liquid and some gas properties.

From the molecular theory of collisions, the temperature is nothing but the average kinetic energy of the molecules in the system, and the pressure is created when these moving molecules collide with the boundary of the system. Therefore, when there are significant intermolecular forces at the microscopic level, the average kinetic energy of the molecules, i.e., the temperature of the system, may increase or decrease, and similarly, the specific volume of the system may be effected. Therefore, a correction factor is needed for temperature and pressure terms in the thermal equations. But any thermal property can be determined by two fundamental quantities from the pressure, temperature, and density. The EOS is the relation that connects these three properties.

Here, $p$ is the pressure (Pa), $ρ$ is the density of mixture (kg/m³), $Rmix$ is the gas constant of that mixture (kJ/kg K), T is the temperature (K), and $Z$ is the compression factor or compressibility factor. For ideal gases, $Z$ is unity and for real gases, it can be either greater or less than one. A $Z$ value less than one represents the attractive forces, and Z greater than one represents the repulsive forces between the molecules. In general, the $Z$ value is less than one after the critical point because as the molecules come closer due to the pressure, they initially experience attractive forces and after a certain combination of pressure and temperature, it again becomes greater than one due to the strong repulsion of electron clouds over the molecules. Also, the behavior of $Z$ varies from species to species, and especially, in combustion it also depends on interspecies interactions. Therefore, the first interesting question to be asked is: Whether sCO₂ operating conditions are in a zone where the molecules are repelling or attracting? This question is answered in the the Compression Factor (Z) section.

In a combustion mixture, there exist polar molecules and nonpolar molecules. The polar molecules are one which will have ionic charges on their outer atoms (for example, H₂O), and the nonpolar molecules will have neutral charges on the outer atoms (for example, CO₂ and O₂). Therefore, to measure the repulsion or attraction forces accurately, one may need to have finer molecular details.

Boyle was the first person to demonstrate the P–T–ρ relation for an ideal gas, but only the extensive work of van der Waals in the late 18th century [18] delineated the first approximation of EOS for real gases. The EOS can be broadly categorized into three types; they are the virial type, molecular-based, and van der Waal type EOS [19,20]. The virial and molecular-based EOS are highly accurate and complex in their formulations. Therefore, using such EOS types in computational fluid dynamics (CFD) combustion applications is difficult because the EOS must consider all the species involved in the mixture and solve for each cell in the computational domain at each time step. The NIST-REFPROP [21] is one such program where complex EOS of such type is used to calculate the thermal properties of the fluids and fluid mixtures. The computational expensiveness involved in using the National Institute of Standards and Technology (NIST) for CFD is reported in Ref. [22]. However, it must be noted that the NIST is considered as the most accurate EOS available. Therefore, it is usual in the literature to see the usage of NIST as a reference for EOS validations where the experimental data are not available. The third category of EOS, i.e., the van der Waals type of EOS, is empirical in nature, and it is the main type used in combustion simulations.

Succeeding proposals after van der Waal have largely modified the basic van der Waals correlation for better accuracies. The improved van der Waals class equations such as Redlich–Kwong (RK) [23], SRK [24], and PRS EOS [25] are the most popular EOS for supercritical CFD simulations of rocket combustion systems [26,27] due to their simple formulation and modest computational cost.

In Eq. (4), the $Z$ can be solved analytically, and it will have three solutions due to its cubic order. Only when the mixture of interest is subcritical in both pressure and temperature, there are three real roots to this equation. Therefore, at subcritical conditions, the correct real $Z$ must be identified by phase equilibrium procedure, and for supercritical conditions, the largest real root can be considered as the compressibility factor.

Here, the subscript $i$ refers to a species index, $Xi$ represent mole fractions, $ai$ and $bi$ are pure species properties, and $kij$ are binary interaction coefficients that are determined empirically. Although empirical in practice, this interaction coefficient is a measure of deviations from the ideal solution behavior for interactions between the ith and jth components [29]. Thus, its value is 1.0 when $i$ equals $j$⁠, i.e., for pure fluid interaction, and it is nearly 1.0 for component pairs which form nearly ideal solutions. Its value differs considerably from 1.0 when the component pair forms highly nonideal solutions. Thus, accurate values of $kij$ are required when $i$ or $j$ is a light hydrocarbon or a nonhydrocarbon (for example, methane with hydrocarbons heavier than n-butane, CO₂–hydrocarbon, H₂S–hydrocarbon, and N₂–hydrocarbon mixtures). More detailed information about $kij$ can be found in Refs. [30–33]. Also, the coefficients $u$⁠, $w$⁠, $a$⁠, and $b$ for Eq. (4) are presented in Table 1. The choice of these coefficients changes Eq. (4) to the EOS of interest.

The CHEMKIN-RG has the capabilities to model the EOSs by van der Waals (vdW), RK, SRK, PRS, Becker–Kistiakowsky– Wilson, and Noble–Abel EOS. In the current work, only the RK, SRK, and PRS EOS are used for comparison, because these are the most popular EOS, which are being used in rocket combustion simulations.

It must be noted that the critical properties such as $Tc$⁠, $Pc$⁠, and $ω$ for all the species and radicals in a combustion phenomenon are not available in the literature. Therefore, it is the practice to assume the critical properties of the largest diluent in the simulation to the species or radicals for which the critical properties are not known.

Since these cubic EOSs are empirical by their origin, adopting them to an application requires a validation with data. For example, some investigations recommend SRK EOS for CH₄/LO_x and kerosene/LO_x mixtures [8,34], whereas Poschner and Pfitzner [35] recommends PRS for H₂/O₂ mixtures. Nonetheless, the most used EOSs are the SRK and PRS. Figures 2(a) and 2(b) show the P–T–ρ correlations of CO₂ and O₂ with various EOS models. Figure 2(a) illustrates a better accuracy for PRS EOS over SRK and RK EOS when compared with NIST, and this accuracy increases as the temperature increases. Interestingly, the accuracy of SRK EOS also increases with temperature, and beyond 1200 K, the SRK and PRS are not distinguishable. For sCO₂, the average deviations of PRS and SRK EOS with NIST are 0.04% and 1.87%, respectively. Further, Fig. 2(b) shows the better accuracy of SRK EOS over PRS and RK. For sO₂, the average deviations of PRS and SRK EOS with NIST are 1.47% and 0.25%, respectively. The percentage deviations of EOS of sCO₂ and sO₂ illustrate that the PRS EOS, which is accurate for sCO₂, is not as accurate as SRK EOS for sO₂. Therefore, the best EOS for combustion mixture depends on the proportions of the CO₂, O₂, and other mixture constituents. It is known that the proportion of these species vary from unburnt to fully burnt condition and at various turbulence levels.

It is shown that the chemical pathways are influenced by the small-scale turbulence (turbulent dissipation rates, N) [13]. The turbulent dissipation rate characterizes the magnitude of molecular mixing (small-scale turbulence) in a combustion process, where the pockets of high strain change in the local chemical reactions which influence chemical kinetic pathways [36]. Therefore, the proportions of mixture constituents vary between the turbulent regimes of different dissipation rates [12]. Hence, this EOS validation is also carried between two turbulent dissipations values such as N = 10,000 and N = 1. Most of the EOS validation in the previous literature is available only for pure species. However, in the current investigation, the EOS validation is carried for combustion mixtures at various RPVs.

For the comparison of EOSs, the determination of mass fractions at various RPVs has been made with the PCMC code as described previously, and CHEMKIN-RG is used for thermal state prediction. The RPV determines the amount of enthalpy released out of reaction from the total available enthalpy of that mixture. When RPV is equal to zero, that is, the mixture is still unburnt, and when RPV is one, the mixture has released its complete enthalpy content, or in other words, the combustion is complete.

Figure 3 is a schematic of the mixture compositions considered for EOS comparison, and the reference inlet mixture is chosen as shown in Table 2. When RPV is zero, that is, at the inlet unburnt mixture, the mole ratio of CH₄ to O₂ is fixed at 0.5 (stoichiometric), and the number of CO₂ moles is 24. This combination of CH₄, O₂, and CO₂ gives the outlet temperature as ∼1500 K. The PCMC first calculates the equilibrium solution of this inlet condition, and then, it solves the PCMC equation (Eq. (1)) between these inlet and equilibrium solutions for various turbulent dissipation values (N). Therefore, based on N value, the proportion of mixture constituents varies between RPV = 0 to 1. The PCMC can calculate the mass fraction of all the species involved in the chemical mechanism, i.e., Aramco 2.0, by a constant enthalpy and constant pressure process. But NIST can only calculate the properties of the mixture having seven species such as CO₂, CH₄, O₂, H₂, H₂O, CO, and C₂H₆. Despite such limited number of species, the NIST can be still used for sCO₂ combustion EOS validation because at any given RPV, the sum of all these mole fractions is more than 99.99%. The operating condition one (OP1) as shown in Table. 2 is considered for the EOS comparison. This operating condition considers the approximate boundary conditions for sCO₂ combustors, i.e., the inlet temperature as 1000 K, exit temperature as 1500 K, and the inlet CH₄–O₂ ratio in stoichiometric proportions.

The difference in the mass fractions of the PCMC solution by SRK and PRS is observed to be negligible (less than 0.001%) for the major species such as H₂O, O₂, and CO₂, but for CH₄, this difference is up to 2.1% and the temperatures differed by 0.12%. Therefore, the average mass fractions and corresponding temperatures of both SRK and PRS solutions have been taken to find the corresponding thermal properties from the NIST. In this comparison, the RK EOS is not considered due to its higher deviation from NIST.

Figure 4 shows the comparison of thermal state prediction of PRS and SRK EOS with NIST at a pressure of 300 atm at two different turbulent dissipation rates. The primary vertical axis represents the mixture density, whereas the secondary vertical axis represents the difference of the EOS prediction with NIST. Due to the variation of species and rise in temperature between the unburnt and completely burnt mixtures at these turbulent dissipation values, a maximum density difference of 2.0% is observed between the real gas EOSs and NIST predictions.

The result shows that the SRK EOS has better accuracy for density prediction at both the turbulent dissipation values. The average deviation of SRK EOS is 0.71%, whereas for PRS, it is 1.78% when N = 1, and at N = 10,000, the SRK EOS is 0.70% and PRS has 1.71%. It must be noted that the uncertainty involved in NIST mixture density calculation is 1%, and SRK EOS solution falls within that uncertainty level. These narrow deviations of SRK EOS with NIST show that SRK EOS is the most accurate EOS which can be used for sCO2 combustors. Also, it is observed that the computational time for these EOS is almost the same because both the EOS are derived from the same cubic equation of state.

When the flow is incompressible, the thermodynamics can be separated from the fluid kinematics (the movement of the flow) and fluid dynamics (the forces in the flow) [37]; it is a simplified assumption for many of the nonreacting flows where the density of the system is almost constant, and this incompressible assumption makes the flow much easier for analyzing. However, the pressure gradients, temperature gradients, and Mach number variations can cause the flow to be compressible. But in supercritical flows, it also appears due to molecular attractive and repulsive forces which are represented by $Z$⁠. Therefore, it is essential to understand the behavior of $Z$ in sCO₂ combustion. In compressible flows, the density variation in the working fluids transfers the energy from the fluid to the surroundings; in other words, the thermal energy in the system converts into mechanical energy (changes in velocity and momentum) or mechanical energy converts into thermal energy (increases the temperature and hence changes in density).

Figure 5 shows the variation in compressibility factor under various operating conditions as listed in Table 2. Here, the operating conditions represent various possibilities of sCO₂ combustor operation. As explained in the Equation of State section, OP1 is the reference case used in the current work, which considers the inlet and outlet boundary conditions for an sCO₂ combustor. The quantification of fundamental thermal properties in this literature is based on the OP1. The OP2 show what may happen to the $Z$ when CO₂ mole fraction increases in the reference mixture, the OP3 shows what may happen to the $Z$ when CH₄ and O₂ mole fractions in reference mixture increases, and the OP4 shows what may happen when the inlet temperature of the reference mixture is decreased. In the later part of this section, an empirical correlation is proposed for the estimation of $Z$ in the sCO₂ combustor. The operating conditions OP5–OP8 are used as additional validating cases for this model.

Also, Fig. 5 answers the question which is being asked in the the Equation of State section, i.e., whether sCO2 operating conditions are in a zone where the molecules are repelling or attracting? The answer is: The sCO₂ combustor exists in an operating zone where there exist repulsive forces among the molecules. The gradual decrement of $Z$ with respect to the progress variable in Fig. 5 shows an important design aspect to the sCO₂ combustor designers. Because, from Eq. (2), the $Z$ is proportional to the pressure, and as the reaction progresses (temperature of the mixture increases), the static pressure may reduce. It may result in the loss of turbine efficiency. In a conventional gas turbine combustor, the loss of static pressure may be due to flow obstacles or due to turbulent mixing; however, in sCO₂ combustors, the designers should deal with the depreciation of $Z$ in the combustor due to the supercritical nature of the flow. It is also interesting to note that, beyond the critical point, up to certain pressures the $Z$ value increases with temperature, and this phenomenon gets reversed after a certain supercritical pressure. The sCO₂ combustor operating conditions are existing in a zone where the $Z$ decreases with temperature.

Figure 5 also shows that the slope of $Z$ is larger for OP3, which indicates that increasing the inlet CH₄ and O₂ would increase the static pressure loss; it is mainly due to the increase in temperature (the Z loss is 5% in this case). On the other hand, OP2 and OP4 show a minimum slope which indicates that the static pressure loss can be minimized either by increasing the CO₂ content in the initial mixture or by decreasing the inlet temperature. It is because the addition of CO₂ absorbs the temperature released due to its high specific heat, and the lower inlet temperature results in a lower final temperature. The dilution of the combustion mixture with additional CO₂ after combustion would help in regaining the $Z$ and hence the static pressure.

In ideal gases, the pressure is nothing but the rate of change of momentum exchange per unit area of the combustor walls. In supercritical conditions, the repulsive forces between the molecules will be added to the overall momentum and hence pressure. Therefore, the pressure correction equation is modified for the pressure implicit split operation algorithm (PISO), and the suitable solution sequence is suggested by Park and Kim [38]. Basically, the $Z$ factor is considered in the PISO algorithm from the real gas EOS, instead of ideal gas EOS.

The calculation of $Z$ is very important because the thermal properties in supercritical combustion are functions of $Z$⁠. Some important isentropic flow relations are presented in Table 3, from the work of Baltadjiev [39]. Therefore, in this work, an empirical model as shown in Eq. (7) is suggested. This equation is derived from the analysis of the operating conditions presented in Table 2. The density values can be computed from the ideal gas assumption and then corrected with the Z calculated from Eq. (7), reducing the computational time. This equation is validated for the inlet pressure 250–300 bar, inlet temperature 800–1200 K, inlet stoichiometric CH₄ and O₂ mixture, and CO₂ mole fraction up to 0.93. Here, $(nCO2)in$ is the inlet CO₂ moles, $(nCH4)in$ is the inlet methane moles, and $Tin$ is the inlet temperature in Kelvin

The calibrated constants are as follows:

$A1=1.0965$⁠, $A2=0.0217,$$A3=0.36e−3$

$A4=0.38e−2$⁠, $A5=0.0272$ , $A6=0.0237$⁠,

$A7=24$⁠, $A8=0.6581$ , $A9=0.00136$⁠, and $A10=1000$

Figure 6 shows the correlation plot between the modeled Z and calculated Z by the SRK EOS. This model is validated with all the operating conditions from OP1 to OP8 in Table 2. Each symbol in the plot represents the Z value at an operating condition over an RPV value from 0 to 1. All these data points in the plot are well inside the 0.5% error lines, which indicates that the proposed model predicts the Z as accurate as 0.5%.

The specific heat of a system indicates the resistance of the system to change its temperature when heat is added. Due to the high heat capacity of CO₂ compared to other species in a combustion mixture, the CO₂ in an sCO₂ combustor carries most the enthalpy with its flow rather than raising the temperature. However, for supercritical combustion, the Z factor also influences the specific heats. An attempt is made in this section to estimate its influence in sCO₂ combustors. Figure 7 shows the variation of $cp$⁠, $cv$⁠, and $γ$ (ratio of specific heats) for combustion mixtures with real gases and ideal gases. The figure illustrates that specific heat capacities calculated for sCO₂ combustion mixtures are always greater than the values calculated by ideal gas assumption because the $Z$ value is always greater than unity. The general formula for specific heat is given by the expression

From Eq. (8), for the same energy content (E) of the system, the specific heats in the sCO₂ combustion system are always higher than the ideal gas case. It must be noted that, as discussed in the the Compression Factor (Z) section, there are some operating conditions beyond the critical point, where the Z is less than unity, and in this region, the specific heats are less than the ideal gas assumption. However, the specific energy in sCO₂ combustion mixture is higher than ideal gas assumption because of the existing repulsive forces. In fact, these repulsive forces when considered increase all the specific properties such as entropy, Gibbs energy, etc.

Figure 8 shows the enthalpy–entropy relation for ideal and real gas assumptions. Here, the enthalpy and entropy are normalized with respect to its initial values at RPV equal to zero. These two curves deviate as they progress in combustion. It explains that there is a higher irreversibility associated with the combustion process due to the influence of Z, and this irreversibility increases at higher temperatures. Also, the enthalpy released is approximately 1.6% more and entropy increased by 0.03% more than the ideal case at the end of the combustion.

The $cp$ of the mixture starts at a value 1.331 and by the end of combustion, it reaches 1.381 kJ/kg K, whereas the $cv$ value starts at 1.111 and reaches 1.187 kJ/kg K. However, it must be noted that the $cp$ and $cv$ values are functions of $Z$⁠, and these values change when the operating conditions change $Z$⁠. Also, the value of $γ$ is about 1.2 for the inlet mixture and 1.164 for the fully burnt products. Another interesting trend can be seen from Fig. 7 that the $cp$ and $cv$ of ideal gas and real gases converges as the reaction progresses, it is because the repulsive forces decrease with increase in temperatures.

As shown in Table 3, the isentropic processes for real gases in the $p−v$ coordinate system are described as $pvns=pvγ/(βTP)=C$ [39]. It is a crucial parameter that determines the path functions such as work and heat involved in a process.

Therefore, in this section, the variation in $ns$ with respect to the reaction progress is discussed. Figure 9 shows that the value of pressure exponent $ns$ is 1.32 at the entry of the combustor, and it gradually reduces to 1.245 by the end of combustion. It is because the $βT$ which is constant for ideal gases is increasing as the combustion progresses, and it results in a reduction of $ns$⁠. In general, for hand calculations, it is very common to assume the process index to be constant, but in sCO₂ applications, this assumption may result in huge errors because the change in $ns$ is up to 7% between the unburnt and fully burnt mixture.

The isothermal compressibility $βT$ is the reciprocal of bulk modulus $(K)$ of the combustion mixture. At first, the $βT$ in the combustion mixture is less than its ideal gas value because the existing repulsive forces between the molecules resist the bulk compression. Second, it increases with respect to the RPV because the repulsive forces decrease as seen in Fig. 5, and hence, the combustion mixture becomes relatively more compressible but still resistive to bulk compression compared to the ideal gas assumption. Therefore, the isothermal compressibility will be overestimated up to 7% when the real gas assumption is not used.

The speed of sound determines how fast a disturbance travels in a flow or, in other words, how fast information of an obstacle spreads in the flow. After a certain ratio of the flow speed and sound speed, the molecules will tend to compression or rarefaction and result in the significant density changes in the flow. Therefore, an attempt is made in this section to estimate the speed of sound in the sCO₂ combustor. As shown in Table 3, the speed of sound is a function of $ns,Z$⁠, and $T$⁠. Therefore, $ns$ and $Z$ determine the speed of sound in the sCO₂ combustor. From the previous discussions, the values of these three parameters are larger when calculated with the ideal gas assumption.

Therefore, the speed of sound must be higher. Figure 10 shows the comparison of sound speed between the real gas assumption and ideal gas assumption. Also, the existing relation between the speed of sound and density fluctuations in a flow has been modified for real gases.

From Table 3, $a2=nsZRT$$⇒dP=a2∗dρ$⁠.

The term $a2$ in Eq. (12) is higher for an sCO₂ combustor than the ideal gas value (⁠$a2=γRT$ in the ideal gas assumption). Therefore, for a velocity change, the change in corresponding density is lower for sCO₂ combustion by a factor $γ/(nsZ)$⁠, which is equal to ∼13%.

Modeling of viscosity is an important aspect in the simulation of any fluid flow where the shear stresses are prominent. Though the viscosity of sCO₂ has huge variation near to critical point, in the operating regime of sCO₂ combustor, as the temperature increases, the viscosity follows the linear trend like an ideal gas, and the viscosity of sCO₂ is higher compared to the viscosity of CO₂ at atmospheric pressures. This section gives an overall view of various viscous models which are being used in current supercritical mixture viscosity formulations and the best fits for 1D and CFD combustion modeling.

Figure 11 shows the flowchart of various methodologies which are being used currently to calculate viscosity in supercritical CFD simulations. There are basically two approaches, shown as two columns in Fig. 11, to estimate the viscosity of the supercritical mixture. One is to calculate the individual species viscosity at required pressure and temperature and use mixing rules to calculate the resultant mixture viscosity [5]. It is represented in the left column of Fig. 11. To calculate the individual species viscosity in the mixture, there are the Lucas method [40], which needs the temperature and pressure as input variable, and Chung et al. method [41], which needs the temperature and density as inputs. Out of these two methods, one can use either the Lucas method [40] or the Chung et al. method [41] for the viscosity calculation of the species involved in the sCO₂ combustion mixture and then use these species viscosities in the mixing rules to find the mixture viscosity. Though the Wilkes mixing method is defined for calculating the viscosity of low-pressure mixtures, it has been used for the supercritical simulations because it is a subset of Wassiljewa–Mason–Saxena (WMS) method [5]. The WMS is a thermal conductivity modeling method. Both the Wilkes method and WMS method have a common variable, and using both together provides consistency. The other column in Fig. 11 shows the one fluid approach. In this approach, the mixture is assumed as one fluid, and all the critical properties of the species are combined to find the effective critical properties, and the viscosity is calculated based on these effective critical properties of the fluid. The Chung and Lucas methods are two popular one-fluid mixture viscosity methods for supercritical applications [42].

A systematic investigation is presented in this section to estimate the suitable viscosity model for sCO₂ applications. Initially, the viscosities of prominent species are compared between Chung, Lucas, and NIST methods. Further, the mixture viscosity methods are compared at various combustion progress variable values. Since there is no reference mixture viscosity to compare the models, the suitable mixture models for sCO₂ combustion are judged based on the computational time, complexity, and performance of the mixture models when applied to individual species. This plot has viscosities of CO₂, H₂O, CH₄, O₂, and CO at 300 bar pressure. Here, the viscosities of CO₂, H₂O, and O₂ are shown between the temperature limits 600–1600 K, and for CH₄ and CO, the temperature is shown only up to 600 K and 500 K, respectively. It is because the NIST data for these species are available only to these temperatures.

Figure 12 shows the prominent species viscosity modeled with Lucas and Chung methods and the comparison of those models with NIST. The CO₂ viscosity plot has an additional viscosity model called Heidaryan et al. [43]. This model is validated between the pressures 75 bar and 1014 bar and between the temperatures 310 and 900 K. It is a simple empirical correlation which only needs pressure and temperature as an input. The plots show that, for CO₂ and O₂, the Chung et al. model is closer to NIST than the Lucas et al. model. The Heidaryan et al. is closer to NIST up to 900 K and deviates later. However, for H₂O, the Lucas et al. model is closer to NIST. Further, within the available NIST data limits, Lucas et al. is closer to NIST for CH₄, and Chung et al. is closer to NIST for CO. Therefore, Fig 12 shows that a best individual species viscosity model varies from species to species. To identify more suitable species viscosity model for sCO₂ combustion applications, the typical weightage of this species distribution is considered, and the deviation from NIST is calculated based on this weighted deviation. Table 4 shows the comparison of these models and their weighted deviation from NIST. The distribution of these species is estimated by using the CMC method as discussed earlier. The average of weighted deviation of Lucas et al. is 0.72%, whereas for Chung et al., it is 0.74%. There is not much difference between both the models. Therefore, the computational time and complexity in formulations determine the suitable viscosity model. It must be noted that CH₄ and CO are not considered for weighted deviation estimation because the NIST data are unavailable for those species at the operating temperatures of sCO₂ combustors.

As shown in Fig 11, the Wilkes method is a low-pressure mixture viscosity method. However, it is also used in supercritical CFD simulations [5] to calculate the mixture viscosity. The more complete information about these models can be obtained from the references provided [42]. There is a simpler model suggested by Brokaw [44] for nonpolar gaseous mixture where the inputs are just the gas composition, viscosities, and the molecular weights of the constituents. It is a derivative of molecular-based Chapman–Enskog theory. This method claims its applicability for high pressures beyond the critical point where there is no large concentration of free radicals. Also, the largest error reported with this mixing rule is less than 4%. Since the major content (more than 90%) of the sCO₂ combustion mixture has nonpolar molecules such as CO₂, CH₄, and O₂, and the expectable radical concentration is comparatively low due to high specific heat of CO₂, this particular model [44] would be an accurate alternative fit for sCO₂ combustion applications, but needs more investigation because the combustion products contain H₂O which is a polar molecule.

where $mij=[4MiMj/(Mi+Mj)2]1/4$⁠.

In the current section, the mixture viscosity models are analyzed. As discussed, the Chung et al. and Lucas et al. are two prominent mixture viscosity methods that are currently being used in supercritical combustion literature. One basic difference between these two mixture models is the method of calculating the combined critical properties. The detailed formulations of these models can be found in Ref. [42]. Unlike the Lucas et al. model, the Chung et al. model considers the interaction of one species with another species while calculating the combined critical properties. Therefore, in the Chung et al. method, if we have “k” number of species in combustion mechanism, there will be a k × k matrix for calculating each combined thermal property of the system. Hence, computationally, the Chung et al. is expensive. Also, the Chung et al. method needs the binary interaction parameter called “k_ij.” As discussed in “The Equation of State” section, the binary interaction parameter describes whether an interaction between two species is ideal or nonideal. This parameter is empirical, and there are many models as described earlier which calculates this parameter. However, these models are not evaluated with experiments for all the species involved in combustion. Therefore, it is usual practice to assume the unknown k_ij as unity. In the current section, the Chung et al. mixture viscosity method is evaluated by considering both k_ij as unity and by choosing appropriate values for it.

Figure 13 compares the viscosity of an sCO₂ combustion mixture calculated by using various mixture models at three different temperatures: 1000, 1250, and 1500 K. These three temperatures are corresponding to unburnt, half-burnt, and fully burnt conditions of the sCO₂ mixture. Again, the species distribution is solved by using the PCMC code at a dissipation rate of unity. Figure 13 also shows the viscosity calculated with Wilkes model and a weighted average. These two models require the individual species viscosities, and they are calculated using the Lucas method of individual species viscosity. At all the temperatures, Lucas et al. has the highest calculated mixture viscosity, and Chung et al. has the lowest. The Chung et al. when k_ij is unity predicts slightly higher viscosity than Chung et al. At 1000 K, a 10% deviation line is drawn below the highest viscosity model, i.e., the Lucas et al. to estimate the deviation of other models. This deviation shows that all the models used here are predicting the mixture viscosity within a 10% deviation. Here, it must be noted that there is no reference data to compare these models to predict the best suitable viscosity mixture model. However, in this section, it is seen that the Lucas et al. is predicting the species viscosity more accurately than other models (when CO and CH₄ predictions are taken into account along with other species). Therefore, though there is no standard reference to identify the superiority of one model with other, the Lucas et al. model can be considered as accurate due to its performance in calculating individual species viscosity. The dynamic viscosity may not be an influencing parameter in a pure turbulent regime; however, when it comes to CFD, the viscosity term in the diffusing term is important. Also, the viscosity is more important in simulating the internal flows, viscous sublayer, Prandtl, and Schmidt numbers. In turbulent flows, the viscosity damps the small scales; therefore, the simulations such as direct numerical simulation (DNS) and large eddy simulation (LES) to capture small scales will be highly effected by the viscosity model. Therefore, a careful judgment must be made in choosing the viscosity model for DNS and LES simulations.

Figure 14 compares the computational time involved in calculating viscosity by these models. These models are simulated in matlab, and the time is multiplied with a factor of 10,000. It is to show the computational time if there are 10,000 cells in the computational domain. The plot shows that Lucas et al. is least expensive in terms of computation, and Chung et al. is more expensive. It is because the Chung et al. model accounts for the interaction of species, and it adds many two-dimensional matrices in the computation. The Chung et al. is much more expensive for LES, DNS, and detailed kinetic simulations. The Wilkes mixing rule and weighted averaging is more expensive than Lucas et al. because they need the input of individual species viscosities. Therefore, Lucas et al. is recommended for sCO₂ simulations due to its lower computational time and performance in calculating the individual species viscosity.

Thermal conductivity is the capacity of a system to transfer the heat energy from it. The thermal conductivity of the system increases as the pressure increases due to the increase in the density. However, it is important to note that, as the temperature increases, at low pressure, the thermal conductivity tends to increase, but beyond the critical pressure and up to certain pressure limit, the thermal conductivity reduces with temperature [42], and after this pressure limit, the thermal conductivity increases again with the temperature.

Figure 15 shows the current standard practices for calculating thermal conductivities in supercritical combustion simulations. The Chung et al. model for calculating thermal conductivity is the most widely used model. The Chung et al. high-pressure thermal conductivity model needs the low-pressure mixture viscosity as an input to estimate the high-pressure thermal conductivity. The first two columns in Fig. 15 illustrate the calculation of thermal conductivity by Chug et al. high-pressure model. The first column shows that the mixture viscosity at 1 atm is estimated by Wilkes methods, and the input of individual species viscosities for Wilkes method is calculated either by the low-pressure Chung et al. or Lucas et al. methods. The second column in Fig. 15 also shows the estimated thermal conductivity by the high-pressure Chung et al. method. However, here the low-pressure mixture viscosity is estimated by one-fluid approach, i.e., either by the Chung et al. method or the Lucas et al. method. The third column shows the proposed method, the estimation of supercritical thermal conductivity by the Stiel and Thodos method. The Stiel and Thodos [42,45] is a simple, generalized correlation for a high-pressure gas thermal conductivity of a mixture. It says that the excess thermal conductivity (excess from low-pressure thermal conductivity) associated with a gas or mixture is a function of its critical properties, density, and molecular weight. This correlation has been validated against 20 datasets of nonpolar gases including CO₂. This formulation is used for calculating the thermal conductivity of species involved in supercritical mixing studies [5]. One input to this Stiel and Thodos model is the low-pressure thermal conductivity of the mixture. The low-pressure thermal conductivity of the mixture can be estimated by the Wassiljewa model (as shown in Eq. (15)). Also, this method is used by Masi et al. [5] for supercritical multicomponent mixing simulation due to its simple formulation, though this model does not take any supercritical condition into account

where x_i and x_j are the mole fractions of the species i and j, and λ_i is the thermal conductivity of the species i.

The bracketed term in $Bij$ also appears in Wilke's method for mixture viscosity correlation. The $ε$ value is suggested as 1.065, later its value is modified by Tandon and Saxena to 0.85 [42].

It must be noted that Masi et al. [5] used Wassiljewa method for calculating the high-pressure thermal conductivity, though it is defined for low-pressure thermal conductivities, and the required individual species viscosities for this model are calculated by the high-pressure Lucas et al. model. However, in the proposed approach, the individual species viscosities are estimated by the low-pressure Lucas et al. model over a wide range of temperatures at 1 atm pressure and tabulated. Later, the tabulated data are called into the Wassiljewa thermal conductivity model to estimate the low-pressure mixture thermal conductivity of the required mixture. Further, the Stiel and Thodos method is used for calculating the supercritical mixture thermal conductivity.

There is no reference mixture thermal conductivity data to identify the superiority of the Stiel and Thodos method over the Chung et al. method. However, when these models are used to identify individual species viscosities, the Stiel and Thodos method shows a better accurate match with NIST data. Figure 16 shows the comparison of individual species viscosities between Chung et al., Stiel and Thodos, and NIST. The Amooey et al. is an empirical thermal conductivity model for CO₂ [46]. Therefore, the Amooey et al. model is also shown in CO₂ thermal conductivity plot in Fig. 16. It shows that the Amooey et al. is accurate only up to 900 K and deviates significantly after that. The Stiel and Thodos method predicts the CO₂ and O₂ thermal conductivities very accurately. For H₂O, it is more accurate than the Chung et al. The predictions of CH₄ and CO are better with Chung et al. than Stiel and Thodos.

Table 5 shows the weighted deviation of the thermal conductivities calculated by the models from the NIST data at 1250 K. The weighted average deviation for Stiel and Thodos is 0.91%, whereas for Chung et al., it is 9.15%. It shows that Stiel and Thodos thermal conductivity model is more accurate for sCO₂ combustion simulations.

Figure 17 shows the comparison of mixture thermal conductivity calculated by Chung et al. and Stiel and Thodos methods. The Chung et al. thermal conductivity model used here is corresponding to the second column of Fig. 15. Both these methods show the increase in thermal conductivity with temperature. Also, the Chung et al. predicts the thermal conductivity better than the Stiel and Thodos method by 15%. The 15% difference in thermal conductivity may affect the heat energy transfer significantly. Besides the better performance of Stiel and Thodos for individual species thermal conductivities, it is also proven to be computationally modest. Figure 18 shows the computational time comparison between both these methods. The Stiel and Thodos method shows the least computational time because all the necessary individual species viscosities are tabulated over the required temperature range. The thermal conductivity simulations are carried out in matlab, and the necessary species distribution with respect to temperature is identified by the PCMC code. Here, the computational times are multiplied with 10,000. The difference between both the approaches will be significant for LES, DNS, and simulations with many species in the chemistry.

Here, the available relations for EOSs and other thermal properties needed to model supercritical CO₂ combustion were reviewed and recommendations were given.

The cubic EOS must be validated for each application due to its empirical nature. The PRS EOS better predicts the thermal state of pure CO₂ compared to SRK and RK. However, over 1000 K, the PRS and SRK are indistinguishable. But, for sCO₂ combustion mixtures, in all the turbulent regimes, the SRK and PRS EOS predict the densities by 0.7% and 1.71% when compared with NIST; which concludes that the SRK is the most appropriate EOS for sCO₂ mixtures.

The projected sCO₂ operating conditions for a supercritical regime where there always exist the repulsive forces (which are quantified by compression factor $Z$⁠) between the molecules and these repulsive forces decreases as the combustion reaction progresses from unburnt to the fully burnt condition. The sCO₂ combustors are expected to observe a new kind of static pressure loss due to the reduction of $Z$⁠. This reduction is more when there is more inlet CH₄ and O₂ mixture, and it can be reduced by decreasing the inlet temperature or increasing the inlet CO₂ moles. The fact that the compressibility factor $Z$ is always greater than unity signifies that the repulsive forces between the molecules always exist, and there exists an energy transfer from the working fluid to the system boundaries.

Also, a new empirical model for $Z$ is proposed for predicting its value at any stage of the reaction progress. This model has been validated with $Z$ calculated with SRK EOS over a wide range of sCO₂ combustor operating conditions. This model can be used with 1D sCO₂ combustion system simulations or in CFD.

Unlike the ideal gases, in real gases, the pressure exponent $ns$ differs from the ratio of specific heats $γ$⁠. The value of $γ$ is about 1.2 for the inlet mixture and 1.164 for the fully burnt products. The resistivity of the sCO₂ mixture to change its temperature is higher due to the existing intermolecular repulsive forces. In fact, all the specific thermal properties are expected to be higher due to the decrease in density. For example, the enthalpy increases by ∼1% and entropy by 0.5%. The pressure exponent changes from 1.32 to 1.24 between the unburnt and fully burnt mixture. Further, the sCO₂ combustion mixture also becomes resistive to the compression by 6.5–9% when compared to the ideal gas assumption. The speed of sound is higher in sCO₂ combustion by a factor $[γ/(nsZ)]1/2$⁠. The combustion mixture is more resistant to density fluctuation induced due to changes in the velocity field by a factor of $γ/(nsZ)$⁠.

The popular viscosity and thermal conductivity models are evaluated based on their accuracy and computational cost. The results show that the Lucas et al. mixture method is more suitable for viscosity modeling, and the Stiel and Thodos method is suitable for thermal conductivity modeling of the sCO2 mixture.

This work was supported by the Department of Energy. This report was prepared as an account of work sponsored by an agency of the U.S. Government. Neither the U.S. Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agency thereof.

• Office of Fossil Energy (Grant No. DE-FE0025260).

a =

mixture a (dyne cm⁴/mole²)

A* =

dimensionless form of “a” for mixture

A1–A10 =

model constants for Z

b =

volume correction term in cubic equations of state (cm³/mole)

b_m =

mixture b (cm³/mole)

B* =

dimensionless form of “b” for mixture

CHEMKIN-RG =

CHEMKIN for real gases [11]

CMC =

conditional moment closure

EOS =

equation of state

f_ω =

function of acentric factor

N =

turbulent dissipation rate

n_s =

process index

NIST =

National Institute of Standards and Technology

N|ζ =

conditioned scalar dissipation

P =

pressure

P_c =

critical pressure (dyne/cm²)

PCMC =

premixed conditional moment closure code [13]

PRS =

Peng–Robinson equation of state

$Q̃i$ =

conditioned mass fraction

R =

universal gas constant (ergs/mole K)

R_mix =

gas constant of the mixture (kJ/kg K)

RK =

Redlich–Kwong equation of state

RPV =

reaction progress variable

S_c =

source term (1/s)

$S̃c|ζ$ =

conditioned source term for RPV

sCO₂ =

supercritical CO₂

sO₂ =

supercritical O₂

SRK =

Soave–Redlich–Kwong equation of state

T =

temperature

u =

a constant in the cubic equations of state

vdW =

van der Waals equation of state

w =

a constant in the cubic equations of state

X_i =

mole fraction of species “i”

Y_i =

mass fraction of species i

Z =

compression factor

$βP$ =

isobaric compressibility

$βT$ =

isothermal compressibility

λ_i =

thermal conductivity of species i

λ_mix =

thermal conductivity of the mixture

μ_i =

viscosity of the species i

μ_mix =

viscosity of the mixture

$ρ¯$ =

density

ω =

acentric factor of the species

$ω̃i˙|ζ$ =

conditioned reaction rate