Estimation of Thermochemical Properties of Propellants Using Peng-Robinson Equation of State

Several aspects make nitrocellulose based powders adequate as a solid propellant for rockets and missiles and the prefe rable propellant to be used in firearms and artillery: it produces less smoke and less fouling than other propellants; its burning rate shows a well-defined relation to pressure and the composition of the gaseous mixture resulting from its decomposition by burning can be predicted with relatively high levels of accuracy. This last feature, if coupled with adequate thermodynamic models, leads to a reliable framework for the simulation of internal ballistics and, ultimately, to an efficient design. Some simpler models use an ideal gas approach, which is obviously not adequate, and some others perform corrections based on the Virial E quation of S tate for sake of algebraic simplicity. This work employs the Peng-Robinson equation of state to model the pressure-volume-temperature (PVT) behavior of the gases produced during the burning reaction, once it is known to be accurate in the high pressures obtained inside combustion chambers. The results were compared to experimental data obtained in a closed vessel device and showed that the Peng-Robinson equation of state could predict the chamber pressure with higher accuracy.


INTRODUCTION
Propellants are a class of explosives that, when properly initiated, can burn in a rapid and controlled manner, a process termed deflagration, thus generating enough gaseous products to propel a weapon projectile (Venugopalan 2015).Despite the drawbacks of producing a large volume of soot and smoke during its combustion, black powder was the first gun propellant ever made, being used for hundreds of years.From the nineteenth century on, this scenario came to be profoundly modified with the discovery of nitrocellulose (NC).This newfound substance became the main ingredient of propellants, resulting in an improved performance, with practically no combustion residue.Due to this feature, these propellants came to be called smokeless powder.Even though they are called powders, they are designed as grains that display different geometric shapes, 2 such as: flakes, ribbons, spheres, cylinders or tubes.Every shape will present a specific burn area which will affect the burn ratio.Along with a lot of possibilities regarding grain shape and size, the smokeless propellants also vary conforming to the ingredients selected.
Conventional propellants are classified according to its composition, mainly the energetic base.When there is only NC as the energetic base, the propellant is called a single-base (SB) type.It can be employed in several systems, from handguns, as pistols and revolvers, to artillery weapons, rockets and missiles.When the propellant is made of NC and nitroglycerine (NG) as the energetic base, it is said to be double-base (DB) type.The presence of a second ingredient makes them more energetic, however a proper care regarding the percentage of NG must be taken, since too much energy makes the gases hotter than necessary, causing erosion to nozzles and barrels, thus reducing its service life.DB applications involve some pistol ammunition, but mostly mortar augmenting charges and military rockets and missiles.Triple-base (TB) propellants receive this classification for being made of NC, NG and nitroguanidine (NGu).NGu is a component whose main contribution is to decrease the adiabatic flame temperature; thus, a TB propellant is applied for large caliber naval guns and artillery weapons with a high rate of fire.Lastly, there is the multibase or nitramine propellant, which is an evolution of the triple base, where NGu is partly or completely replaced by RDX (cyclotrimethylenetrinitroamine).
In terms of energy, a nitramine propellant remains between TB and DB.
When developing new propellants formulations or designing new improved-performance systems, a main concern must be a more accurate and more reliable estimation of thermochemical parameters, such as projectile muzzle velocity and chamber pressure.Usually, the only input available for the calculation of these parameters are the loading density and the propellant's composition.Therefore, obtaining precise data to implement into internal ballistics' models is a major concern.In order to address this problem, one should take into account mechanical, kinetic, geometric and combustion thermodynamics aspects (Suceska, 1999;H. M. Stationary Office, 1987;Harvazinski and Talley, 2018;Dyer et al., 2007); Congiunti and Bruno, 2003;Hickey and Ihme , 2013) In this paper, a case study of a propellant combustion under closed volume conditions is presented to determine thermodynamic factors that affect internal ballistic.These factors comprehend the composition of combustion products assuming water-gas equilibrium and the calculation of energy balance, adiabatic flame temperature and chamber pressure.
Moreover, a special highlight to the difficulties involving the selection of a proper mathematical model to describe the chemical equilibrium in combustion products must be considered.Many studies approaching this problem have chosen the classical Virial-2 equation of state (EoS) to represent real gases (Saurel and Neron, 2021;Suceska, 1999;Xu et. al., 2016) and, to our best knowledge, this choice is based only on the algebraic convenience (Longdon 1987).However, the classical Virial-2 equation of state presents relevant deviations when modeling high pressures and temperatures phenomena, due to its quadratic form (Prausnitz et al. 2000).To tackle these limitations, a Peng-Robinson Equation of State (PR-EoS), which has a cubic form, is selected, because of its known reliability when modeling high pressure gaseous systems (Bergan 1991;Congiunti and Bruno 2003;Giorgi et al. 2014;Hickey and Ihme 2013;Jofre and Urzay 2021;Mallepally et al. 2019;Müller et al. 2017;Müller and Pfitzner 2017;Myint et al. 2016).Additionally, we report the results of a case study where computer simulations were conducted to predict the combustion chamber pressure using the PR-EoS for a DB propellant formulation.These computer results using PR-EoS were compared to the classical Virial-2 EoS approach and to experimental data obtained from a (ballistic) closed vessel.

METHODOLOGY Experimental Protocol
Even though parameters such as chamber pressure and others are necessary for the design of guns and propellant charges, a gun barrel is not a suitable apparatus for determining transient high pressures and how they change over time due to the rapid forward movement of the projectile, the gas and the burning propellant.For measuring transient pressures, a constant volume condition is mandatory, and the suitable apparatus for that is the closed vessel (Grivell 1982;Mukhtar et al. 2019).
A Closed Vessel Test consisted on enclosing and burning a propellant charge of formulation A within a vessel of constant volume to estimate their ballistic parameters (Oliveira et al. 2005).The equipment is represented in Fig. 1   This test was conducted for six different loading densities (Mehta et al. 2015).After closing the vessel, the charge was ignited by the passage of electric current through a filament that dazzled and ignited a small amount of priming (black powder).
Although black powder was used, it is also very common to use squibs for this purpose.This ignition was remotely controlled by a computer that acquired pressure versus time data.The closed vessel has a chamber volume of 200 cm 3 and the tested charges had loading densities of 0.1, 0.12, 0.14, 0.16, 0.18 and 0.20 g/cm 3 .The charge dimensions were 17.5 cm x 15.3 cm x 0.038 cm and the composition of Formulation A is described in Table 1.Before plotting the pressure versus time data, the equipment manufacturer was contacted to provide the storage structure of the data in the output (binary) file.The encoding provided is shown in Table 2.
With the format provided, a program was built for the decoding of this information and its availability for later statistical treatment was stored in a vector variable termed as P exp , with a component for each of the 3072 time points acquired.

Propellant formulation
The propellant deflagration is traditionally modeled as a decomposition reaction (Cooper and Kurowski 1997;Longdon 1987), which will be described as follows.The main component of DB propellants is NC (Longdon 1987), which can be obtained by the nitration of cellulose as described in Fig. 2.
As it can be seen, the monomer formula can be expressed as C 6 H 7 O 2 (OH) 3-x (ONO 2 ) x , where x can assume values from 1 to 3, which leads to a nitrogen content from 6.76% to 14.14%.Nitrocellulose can be gelatinized by an ether-alcohol mixture, acetone and even NG to generate a material that can be conformed into different shapes.Several substances are also incorporated, such as stabilizers (diphenylamine, ethyl centralite), plasticizers (dibutyl phthalate), flame suppressors (potassium salts) and burn rate modifiers (lead and copper salts), or added to the propellant surface to provide characteristics such as the ability to dissipate static electricity or to act as a solid lubricant (graphite) (Cooper and Kurowski 1997).In general, the substances used on propellants only exhibit carbon, oxygen, hydrogen and nitrogen in their molecules, as can be seen on Table 3 (Rumble 2021).Therefore, it is a common practice to express the reaction mixture as a hypothetical substance, with a formula given by C a H b N c O d , which is obtained as a mass-average of its constituents.

Assumptions
To establish a model for propellant formulation designs, the following assumptions are made (Sudarsan et al. 2016): • Burning of propellant is a continuous process, self-sustained, smooth and layered; • Propellant combustion undergoes breaking of chemical bonds described by elementary reactions in accordance with Kisliakowsky-Wilson (K-W) rules and its modified version; and • Equilibrium constant of water-gas equilibrium is temperature dependent only.

Mathematical model
The following model is developed in three steps: • Step 1 presents the oxygen balance for the explosive species; • Step 2 calculates energy balance and adiabatic flame temperature; • Step 3 gives thermochemical properties of propellant assuming water gas equilibrium and PR-EoS approach; and • Step 4 implements thermodynamics data into a computer simulation framework.

Oxygen balance
If the total combustion is achieved, as represented by the chemical reaction (Eq.1): (1) the oxygen balance (Ω), which is the ratio between the remaining mass of oxygen and the original mass, is given by Eq. 2: (2) It can be shown, from data on Table 3, that NC (-73.4% < Ω < -24.2%) and NGu (-30.7%)exhibit strongly negative values for Ω and only NG (3.5%) exhibits a slightly positive one.In fact, for most propellants, -50% < Ω < -30%, which leads to incomplete combustions, with the generation of H 2 , CO and even C along with the gases shown in Eq. 1.

Energy balance and the adiabatic flame temperature
The standard heat of combustion of propellants (Qv) is usually measured in an adiabatic calorimeter at constant volume, which leads to Eq. 3 (Van Ness et al. 2017): (3) where Σ p ΔU o f is the energy of formation of the gaseous products and Σ r ΔU o f is the energy of formation of the reactants, both evaluated at standard conditions -usually 298 K -and given in Tables 4 and 5 (Longdon 1987).Based on Table 4 and 5, it is possible to estimate the adiabatic flame temperature (T o ), which is the one achieved when the heat of combustion is used to heat the gaseous products, as in Fig. 3.   Assuming that the energy is not a function of pressure, it is straightforward to calculate T o (in K) (Eq.4) based on the isochoric heat capacity (C v ) (in J/mol.k), which can be obtained by Eq. 5 (Van Ness et al. 2017): where R is the universal constant of gases and C p (in J/mol.K) is the isobaric heat capacity.This parameter can be calculated by applying a correlation for the gaseous products (Eq.6) based only on the temperature (in K), where the required coefficients for each gas are given in  Considering that the energy balance and the adiabatic flame temperature depend on the composition of the gaseous mixture obtained, this topic will be addressed in the following section.

Gaseous products and the water-gas equilibrium
It is assumed, at this point, that the combustion product is a mixture of CO, CO 2 , H 2 , H 2 O and N 2 , with the first four being involved in the so-called water-gas equilibrium reaction (Eq.7) (Rumble 2021): (7) whose equilibrium constant is given by Eq. 8 (Smith et al. 2010): (8) where ΔG O is the Gibbs free energy, ΔH O is the enthalpy and ΔS O is the entropy, all at standard condition.Some works on ballistics relate this constant directly to the concentration of the gases involved (CO 2 , H 2 , CO and H 2 O), assuming ideal gas behavior (Prausnitz et al. 2000).Afterwards, a correction is conducted using the Virial equation of state (Van Ness et al. 2017), truncated at its quadratic term (V-2 EoS) (Longdon 1987).The procedure built this way is somehow convenient, since it results in an iterative scheme for T o , with the composition being calculated analytically in each iteration (Prausnitz et al. 2000).Besides, it is algebraically independent of the pressure in the vessel, which is calculated after convergence is achieved (Longdon 1987).
In this work, we opted to use the Peng and Robinson Equation of State (PR-EoS) (Van Ness et al. 2017), illustrated in Eq. 9, which shows to be more accurate for the high pressures found in ballistic cycles: (9) where v is the molar volume, R is the universal constant of gases, T is the temperature and a (T) and b are parameters calculated by mixing rules given in Eqs. 10 and 11, as follows (Prausnitz et al. 2000): (10) (11) where x i is the molar fraction of "i".The parameters for individual species are related to their critical coordinates (T c i and P c i ) and their acentric factor (ω i ) as follows in Eqs. 12 and 13 (Prausnitz et al. 2000): The critical coordinates of gaseous species can be found in Table 7, along with their acentric factors (Prausnitz et al. 2000).
Table 7. Critical coordinates and acentric factors of the products.The equilibrium constant is related to the activities (â i ) of the gases involved in the equilibrium, according to Eq. 14 (Van Ness et al. 2017): (14) where fi is the fugacity of "i" in the mixture, φi is the fugacity coefficient of "i" in the mixture and f o i is the fugacity of "i" in the reference state, which can be considered the fugacity of the ideal gas in the same temperature, leading to f o i = P o for all components, thus leading to Eq. 15. ( The fugacity coefficient can be calculated with Eq. 16 (Van Ness et al. 2017): (16 where V is the volume, R is the universal constant of gases, T is the temperature and Z is the compressibility factor, given by Eq. 17. Using Eqs.9-11 in Eq. 16, one can find the fugacity coefficients for PR-EoS, based on Eqs.18 and 19 (Van Ness et al. 2017). ( The reader must be aware that, even though N 2 is not involved in the equilibrium, it must be considered for the calculation of all fugacity coefficients. Finally, the dependence of the energy to the molar volume (ultimately, to pressure itself) is incorporated, in order to correct Eq. 4. For that, the first law of thermodynamics is used, along with one of Maxwell relations (Van Ness et al. 2017), to give Eq.20: (20 where U 0 is the energy of an ideal gas in the same conditions (emulated for v → ∞), which is a function only of the temperature.From Eq. 9: (21) From equation 10: From equation 12: (23) A detailed deduction of equations from Eqs. 18 to 23 can be found in the Appendix.
The complete framework is ready for numerical implementation at this point, but it is possible to reduce the number of variables (which is convenient for computational efficiency), by expressing the number of moles of each gaseous species in terms of the number of moles of CO 2 , which can be achieved by an atom balance (Eq.24). (24)

RESULTS AND DISCUSSION
When designing new propellants formulations, two important parameters must be investigated: the maximum pressure and the flame temperature.They can reflect the burning characteristics of propellants, because their performance parameters, such as burn rate, vivacity, force constant, co-volume and specific impulse are based on the knowledge of the chemical composition and these two parameters.Thus, the maximum pressure and the flame temperature are important input parameters in numerical simulations to predict performance properties.Therefore, looking into different approaches to improve accuracy in theoretical analysis is of great importance for developing new propellants formulations.
Therefore, this work used PR-EoS to estimate the maximum pressure of the propellant formulation A. The results were compared to the experimental values and the calculated value using classical Virial-2 EoS, in order to evaluate the accuracy of the PR-EoS in the calculations of this thermodynamic parameter.

Calculation Results of Maximum Pressure
In order to explain why the pressure is important for propellants, it is worth briefly mentioning part of the ballistic cycle.When a propellant is ignited, within a few milliseconds, hot gases are produced from the burning surface of each grain, and the pressure in the chamber, that started to rise due to the ignition of the primer, keeps rising rapidly.Since most projectiles are resistant to propulsion, an extra force is required to separate the projectile from the case.This resistance makes the chamber pressure get even higher until it is sufficient to engrave the driving band.After that, the projectile starts to be released.Although the available volume increases as the projectile moves along the bore, since the burning rate is dependent on the chamber pressure, the rate of gas production is faster than the projectile speed inside the bore, so the pressure continues to increase up to its peak, until the projectile has moved some distance.
After that, the pressure drops and by the time the projectile has reached the halfway point of the bore, the propellant is all-burned.
By knowing how this cycle works, it is possible to understand how the pressure dictates whether a shot will occur or not, whether the gases produced will be able to produce the desired thrust, leading to an expected muzzle velocity, and whether there will be damage to the barrel.
In this paper, we calculated the maximum pressure of formulation A under different loading densities.The calculated and reference values are shown in Table 8.From the Closed Vessel Testing, a pressure-time profile for six different loading densities, illustrated in Fig. 4, was obtained for formulation A, where the theoretical values calculated using PR-EoS and V2-EoS were represented in the graph as squares and circles, respectively.The experimental load density is calculated by simply dividing the propellant mass employed by the chamber nominal volume and the experimental maximum pressure is the largest of the 3072 recorded measurements.
By analyzing the data in Table 8, it can be seen that PR-EoS gave more accurate results than V2-EoS for all the six loading densities.Moreover, except for Δ = 0.18 and 0.20 g/cm 3 , the remaining pressures were slightly higher than the experimental ones.On the other hand, all pressures calculated using V2-EoS were lower than the experimental values.
Considering that the main purpose of doing theoretical estimations, within the energetic materials field, is to allow empirical researchers to focus their resources on promising new energetic candidates, which offers good performance and safety, thus avoiding unnecessary costs and possible accidents.Therefore, in terms of pressure, a theoretical result that is lower than the experimental value would be undesirable, because, in the case of Δ = 0.16, one might erroneously select a case material designed to support up to 200 MPa without knowing that the pressure could go up to 214.4 MPa, thus causing an accident.Therefore, for predicting maximum pressure for formulation A, PR-EoS would definitely be the better choice, since its deviations are between 0.1-5.3%.

Calculation Results of the Adiabatic Flame Temperature and the Chemical Composition of Combustion Products
The Adiabatic Flame temperature is a thermodynamic parameter that represents the maximum temperature that combustion products can reach when this temperature was raised under adiabatic conditions as a result of the evolution of heat from the reaction itself (Munir and Anselmi-Tamburini 1989).Since this parameter reflects the energy level of propellants, its assessment contributes greatly to understand the reaction energy release process.
The theoretical approach to estimate the adiabatic flame temperature is based on the selection of a suitable PVT EoS for estimating combustion products.Since PR-EoS has proven to be more accurate than V2-EoS regarding the experimental values of maximum chamber pressure for formulation A, this paper assumed this model would also be suitable for predicting the chemical equilibrium of combustion products.Based on this assumption, both the equilibrium chemical composition of the combustion products and the adiabatic flame temperature were calculated using the PR-EoS model.Table 9 shows the predicted adiabatic flame temperature as well as the predicted molar fractions of the gases in the product of each experiment.From Table 9 and Fig. 5, it was possible to observe how the molar fractions from the combustion products changed when the loading density increased.Since N 2 (g) is not part of the water-gas equilibrium, its molar fraction did not suffer any change.
However, for the remaining gases, it could be noticed, in Fig. 5, that whereas x CO 2 and x H 2 increased when the loading density increased, the opposite happened for x CO and x H 2 O , which means that an increase in the loading density favors the production of CO 2 (g) and H 2 (g).By looking into data from Table 9, another correlation can be observed in Fig. 6.In this case, it was straightforward to see that the temperature increased almost linearly with the increase of the loading density of the propellant, displaying a correlation value of 0.9983.This observation is coherent, because the higher the propellant mass within the same container, the higher will be the gas production, thus leading to higher pressures and temperatures.

CONCLUSION
In this paper, we investigated whether the Peng-Robinson Equation of State (PR-EoS) could predict the maximum chamber pressure of a propellant formulation more accurately than the quadratic virial equation of state (V2-EoS).To assess that, we compared the theoretical values from both models to the experimental results obtained from a ballistic closed vessel test for a propellant formulation under six different loading densities.
The comparison showed that PR-EoS could give more accurate results, since its maximum deviation was 5.3%, whereas for the V2-EoS, it was 14.9%.Besides, since all the pressures predicted by V2-EoS were lower than the measured values, this could lead to a poor selection of case material, thus resulting in an accident.
After proving that the PR-EoS was more suitable for predicting the maximum chamber pressure, this model, based on the water-gas equilibrium, was used to describe the thermodynamic state of gaseous combustion products and then the adiabatic flame temperature was calculated.
From the predicted values, we noticed that a higher loading density favors the production of CO 2 (g) and H 2 (g) and increases almost linearly the adiabatic flame temperature.Both the adiabatic flame temperature and the composition of the gaseous products are crucial to the design of tube/barrel guns or missiles/rockets that use such propellant. : Source: Elaborated by the authors.
by the authors.
by the authors.

Figure 5 .
Figure 5. Behavior of molar fractions in function of the loading density of propellant.

Figure 6 .
Figure 6.Behavior of adiabatic flame temperature in function of the loading density of propellant.

Table 1 .
Chemical composition of propellant formulation A.

Table 2 .
Data structure in PE7 files.
Source: Elaborated by the authors.

Table 3 .
Moles of atoms in 1 gram of substance.

Table 4 .
Energy of formation of reactants.

Table 5 .
Energy of formation of (possible) products.

Table 6 .
Parameters for the C p correlation.

Table 8 .
Maximum pressure on the experiment and predicted chamber pressure as functions of the load density.

Table 9 .
Adiabatic flame temperature and composition of the gaseous product as functions of the load density (Δ).