Project 1 - Parsing NASA thermodynamic data

Aim: To Parse NASA Thermodynamic Data in Matlab and obtain plots of Standard State Enthalpy, Entropy and Specific Heat at Constant Pressure for all the species present in the data file and also calculate molecular weight of all the species.

Historical Background: In 1976 NASA came up with a Computer Program for calculation of Complex Chemical Compositions of species and subspecies involved during combustion [1] by Gordon and McBride. This was freely available and marked the first step in studying the products of combustion involved at different stages of the process. This computer program helped in the progress of Rocket Science at that time. In 1990, Sandia National Laboratories operated for the United States Department of Energy improvised the existing computer program by incorporating over 1000 species and subspecies of combustion. This program was called The Chemkin Thermodynamic Data Base [2] and to this date it is the source code of much of the commercial and opensource software for Combustion Simulations. This computer program is currently available commercially as ANSYS CHEMKIN-PRO. Thus the data we are working on is like foundation or corner stone of any combustion research and it is imperative to understand how this data was formulated.

The CHEMKIN Thermodynamic Data Base: It is a general purpose Chemical Kinetics package (computer code) that is used as a tool for solving Complex Chemical Kinetics problems in reacting flows found in Combustion research. This computer program contains curve-fitted data similar to the one used in NASA Complex Chemical Equilibrium Program [1]. Developed for Combustion Research for small molecule Hydrocarbons that leads to NOx pollutants. The whole purpose is to obtain thermodynamic properties of species and subspecies of combustion process such as Standard State Enthalpy, Entropy of various species involved in Chemically Reacting flows as a function of Temperature only. Thus the whole program assumes the species to be functioning as ideal gases. This is because the temperature range we are working on is far above the critical temperature and we know know if Reduced Temperature Tr = T/Tc > 2 then compressibility factor Z == 1, that is the gas behaves as an Ideal gas.

CHEMKIN is Structured into 4 parts:

1. The Interpreter: First the user writes the symbolic description of the reacting Mechanism. The interpreter reads the reaction Mechanism and then extracts the thermodynamic information of species involved from the Thermodynamic Data base. By specific inputs to the interpreter the user can add/ modify the information in the data base.

2. The Thermodynamic Data Base: The vast majority of the data comes from JANAF Polynomials [3] and for specific Hydrocarbons that are not found here the Date is taken from the Third Millennium Ideal Gas and Condensed Phase Thermochemical Data base for Combustion [4] and finally from Quantum Chemistry Computations.

3. The Linking file: The output is printed out and a linking file is generated that contains all the information on species, reactions involved and mechanism of reactions.

4. Gas-Phase subroutine Library: When the linking file is generated the user is set to use it. It allows the user to return information on elements, species, equation of state, thermodynamic properties, chemical production rates & etc.

Methodology used for the Calculation of Thermodynamic DATA:

Initially high-quality ab initio electronic structure calculations were used but due to the availability of supercomputers advanced techniques like BAC-MP4 Method were incorporated. It combines electronic structure calculations 4th-Order Moller-Plasset perturbations theory (MP4) with Empirically derived Bond Additivity Correction (BAC) factors to account for systematic errors in ab initio calculations.
The curve fitting process uses Constrained Linear Least Square fitting procedures. For every species there are 14 coefficients spanning to two temperature ranges. For temperature <= 1000 deg_C for most species is the middle temperature are low temperature coefficients. For temperature beyond this are high temperature coefficients.

Finally the three thermodynamic properties are obtained as follows:

Cp/R = a1 + a₂T + a₃T² + a₄T³ + a₅T⁴
H/RT = a₁ + a₂/2T + a₃/3T² + a₄/4T³ + a₅/5T⁴ + a₆/T
S/R = a₁ln|T| + a₂T + a₃/2T² + a₄/3T³ + a₅/4T⁴ + a7

File Parsing In Matlab:

File parsing in MATLAB involves reading and interpreting the structured text data from the thermodynamic data file. The file is read line-by-line to extract relevant information, such as species names and their corresponding NASA polynomial coefficients. The parsing process includes identifying and extracting numerical values embedded within the text, handling multiple lines for each species, and managing the temperature ranges for which different polynomial coefficients are applicable. MATLAB’s string manipulation and file I/O functions are utilized to systematically retrieve and process the data. Once the data is parsed, it is used to calculate thermodynamic properties, which are then visualized through plots. This systematic approach allows for efficient handling and analysis of large datasets, providing valuable insights into the thermodynamic properties of various species.



Approach used
The approach taken in this project involves parsing the thermodynamic data from a text file in MATLAB to calculate key thermodynamic properties for various species. The properties of interest include specific heat at constant pressure (Cp), enthalpy (H), and entropy (S). The given NASA polynomial coefficients for each species are used to calculate these properties as functions of temperature. The temperature range is divided into high and low regimes, with specific polynomial coefficients applicable to each range. The calculated properties are then plotted against temperature, and the molecular weight of each species is determined based on its chemical formula. This approach ensures a comprehensive understanding of the thermodynamic behavior of different species over a wide temperature range, which is crucial for applications such as combustion research.



Code with Explanation

1. Specific Heat Function

% Function to calculate specific heat

function [specificHeat] = calcSpecificHeat(coeffHigh1,coeffHigh2,coeffHigh3,coeffHigh4,coeffHigh5,coeffLow1,coeffLow2,coeffLow3,coeffLow4,coeffLow5,temp,midTemp,univGasConst)

% Determine temperature range

% The first 7 coefficients are for high temperatures, the second 7 are for low temperatures

if(temp <= midTemp)

specificHeat = univGasConst .* (coeffLow1 + coeffLow2 .* temp + coeffLow3 .* temp.^2 + coeffLow4 .* temp.^3 + coeffLow5 .* temp.^4);

else

specificHeat = univGasConst .* (coeffHigh1 + coeffHigh2 .* temp + coeffHigh3 .* temp.^2 + coeffHigh4 .* temp.^3 + coeffHigh5 .* temp.^4);

end

end

Code Explanation
- This function calculates the specific heat at constant pressure (Cp) for a given temperature.
- It uses different sets of polynomial coefficients for high and low-temperature ranges.
- For temperatures less than or equal to the global mid-temperature , the low-temperature coefficients are used.
- For temperatures greater than , the high-temperature coefficients are applied.

2. Enthalpy Function

% Function to calculate enthalpy

function [enthalpy] = calcEnthalpy(coeffHigh1,coeffHigh2,coeffHigh3,coeffHigh4,coeffHigh5,coeffHigh6,coeffLow1,coeffLow2,coeffLow3,coeffLow4,coeffLow5,coeffLow6,temp,midTemp,univGasConst)

% Determine temperature range

% The first 7 coefficients are for high temperatures, the second 7 are for low temperatures

if(temp <= midTemp)

enthalpy = (univGasConst*temp).*((coeffLow1 + (coeffLow2.*temp)/2 + (coeffLow3.*(temp).^2)/3 + (coeffLow4.*(temp).^3)/4 + (coeffLow5.*(temp).^4)/5 +(coeffLow6./temp)));

else

enthalpy = (univGasConst*temp).*((coeffHigh1 + (coeffHigh2.*temp)/2 + (coeffHigh3.*(temp).^2)/3 + (coeffHigh4.*(temp).^3)/4 + (coeffHigh5.*(temp).^4)/5 + (coeffHigh6./temp)));

end

end

Code Explanation

- This function calculates the enthalpy (H) for a given temperature.
- It distinguishes between high and low-temperature ranges using different polynomial coefficients.
- The enthalpy is calculated based on the temperature and the corresponding coefficients for each range.

3. Entropy Function

% Function to calculate entropy

function [entropy] = calcEntropy(coeffHigh1,coeffHigh2,coeffHigh3,coeffHigh4,coeffHigh5,coeffHigh7,coeffLow1,coeffLow2,coeffLow3,coeffLow4,coeffLow5,coeffLow14,temp,midTemp,univGasConst)

% Determine temperature range

% The first 7 coefficients are for high temperatures, the second 7 are for low temperatures

if(temp <= midTemp)

entropy = univGasConst .* (coeffLow1 .* log(temp) + coeffLow2 .* temp + coeffLow3 .* (temp.^2) / 2 + coeffLow4 .* (temp.^3) / 3 + coeffLow5 .* (temp.^4) / 4 + coeffLow14);

else

entropy = univGasConst .* (coeffHigh1 .* log(temp) + coeffHigh2 .* temp + coeffHigh3 .* (temp.^2) / 2 + coeffHigh4 .* (temp.^3) / 3 + coeffHigh5 .* (temp.^4) / 4 + coeffHigh7);

end

end

Code Explanation
- This function calculates the entropy (S) for a given temperature.
- It uses different sets of coefficients for high and low-temperature ranges.
- Entropy is calculated using the temperature and the appropriate coefficients for the respective range.

4. Molecular Weight Function

% Function to calculate molecular weight

function molWeight = calcMolWeight(speciesName)

% Store species name

name = (speciesName);

% Create an array of elements

elements = ['C','H','O','N','A'];

% Store the atomic masses of the elements

atomicMasses = [12.0107 1.00794 15.994 14.0067 39.948];

% Initialize molecular weight to 0

molWeight = 0;

% Loop through the species name to calculate molecular weight

for i = 1:length(name)

for j = 1:length(elements)

% Compare each character of the species name with the elements array

if strcmp(name(i),elements(j))

% Calculate the molecular weight if there is a match

molWeight = molWeight + atomicMasses(j);

% Store the position for future use

pos = j;

end

end

% Convert string to number

k = str2double(name(i));

% If k is greater than 1, calculate molecular weight accordingly

if(k > 1)

molWeight = molWeight + atomicMasses(pos) * (k - 1);

end

end

% Print the molecular weight of the species

fprintf('Molecular Weight of the species %s is %f\n',name, molWeight)

end

Code Explanation
- This function calculates the molecular weight of a species based on its chemical formula.
- It iterates through the species name, identifying elements and their quantities, and sums up their atomic weights to get the molecular weight.

5. Main Program

% Clear previous data, close figures, and clear command window

clear all

close all

clc

% Open and read the file

fileID = fopen('THERMO.dat','r');

firstLine = fgetl(fileID);

% Define the universal gas constant

univGasConst = 8.314;

secondLine = fgetl(fileID);

% Read the header for global temperature data

B = strsplit(secondLine, ' ');

% Parse global temperature range from the header

globalLowTemp = str2double(B{2});

globalMidTemp = str2double(B{3});

globalHighTemp = str2double(B{4});

% Define temperature range

temperatureRange = linspace(globalLowTemp, globalHighTemp, 1000);

% Skip the next three comment lines

for j = 1:3

tempLine = fgetl(fileID);

end

% Initialize molecular weight array for 53 species

molWeightArray = zeros(53, 1);

% Loop through each species to read and process their data

for i = 1:53

% Parse the species block

line = fgetl(fileID);

splitLine = strsplit(line, ' ');

species = splitLine{1};

% Extract coefficients from the file

line2 = fgetl(fileID);

coeffIndex = strfind(line2, 'E');

coeffHigh1 = str2double(line2(1:coeffIndex(1)+3));

coeffHigh2 = str2double(line2(coeffIndex(1)+4:coeffIndex(2)+3));

coeffHigh3 = str2double(line2(coeffIndex(2)+4:coeffIndex(3)+3));

coeffHigh4 = str2double(line2(coeffIndex(3)+4:coeffIndex(4)+3));

coeffHigh5 = str2double(line2(coeffIndex(4)+4:coeffIndex(5)+3));

line3 = fgetl(fileID);

coeffIndex = strfind(line3, 'E');

coeffHigh6 = str2double(line3(1:coeffIndex(1)+3));

coeffHigh7 = str2double(line3(coeffIndex(1)+4:coeffIndex(2)+3));

coeffLow1 = str2double(line3(coeffIndex(2)+4:coeffIndex(3)+3));

coeffLow2 = str2double(line3(coeffIndex(3)+4:coeffIndex(4)+3));

coeffLow3 = str2double(line3(coeffIndex(4)+4:coeffIndex(5)+3));

line4 = fgetl(fileID);

coeffIndex = strfind(line4, 'E');

coeffLow4 = str2double(line4(1:coeffIndex(1)+3));

coeffLow5 = str2double(line4(coeffIndex(1)+4:coeffIndex(2)+3));

coeffLow6 = str2double(line4(coeffIndex(2)+4:coeffIndex(3)+3));

coeffLow14 = str2double(line4(coeffIndex(3)+4:coeffIndex(4)+3));

% Calculate thermodynamic properties

entropy = calcEntropy(coeffHigh1,coeffHigh2,coeffHigh3,coeffHigh4,coeffHigh5,coeffHigh7,coeffLow1,coeffLow2,coeffLow3,coeffLow4,coeffLow5,coeffLow14,temperatureRange,globalMidTemp,univGasConst);

specificHeat = calcSpecificHeat(coeffHigh1,coeffHigh2,coeffHigh3,coeffHigh4,coeffHigh5,coeffLow1,coeffLow2,coeffLow3,coeffLow4,coeffLow5,temperatureRange,globalMidTemp,univGasConst);

enthalpy = calcEnthalpy(coeffHigh1,coeffHigh2,coeffHigh3,coeffHigh4,coeffHigh5,coeffHigh6,coeffLow1,coeffLow2,coeffLow3,coeffLow4,coeffLow5,coeffLow6,temperatureRange,globalMidTemp,univGasConst);

% Calculate molecular weight

molWeightArray(i) = calcMolWeight(species);

% Save molecular weight to a text file

molWeightFile = fopen('Molecular_Weight.txt', 'w');

fprintf(molWeightFile, 'Molecular Weight of the species %s is %f\n', species, molWeightArray(i));

fclose(molWeightFile);

% Create a directory for the species and save plots

currentDir = pwd;

mkdir(species);

cd(species);

% Plot enthalpy vs temperature

figure(1)

plot(temperatureRange, enthalpy, 'linewidth', 3, 'color', 'b');

xlabel('Temperature [K]');

ylabel('Enthalpy [KJ/mol]');

grid on;

title(sprintf('Enthalpy vs Temperature for %s', species));

saveas(figure(1), 'Enthalpy.jpg');

% Plot entropy vs temperature

figure(2)

plot(temperatureRange, entropy, 'linewidth', 3, 'color', 'r');

xlabel('Temperature [K]');

ylabel('Entropy [KJ/mol-K]');

grid on;

title(sprintf('Entropy vs Temperature for %s', species));

saveas(figure(2), 'Entropy.jpg');

% Plot specific heat vs temperature

figure(3)

plot(temperatureRange, specificHeat, 'linewidth', 3, 'color', 'g');

xlabel('Temperature [K]');

ylabel('Specific Heat [KJ]');

grid on;

title(sprintf('Specific Heat vs Temperature for %s', species));

saveas(figure(3), 'SpecificHeat.jpg');

% Return to the original directory

cd(currentDir);

end

Coding Logic

1. File Reading and Initialization:
- The code begins by opening the thermodynamic data file and reading the first line, which typically contains header information.
- The global temperature range (low, mid, and high temperatures) is extracted from the second line.

2. Skipping Comment Lines:
- The next three lines, which are comments, are skipped using a loop.

3. Main Loop for Parsing and Calculations:
- The main loop runs 53 times (one for each species in the file).
- For each species:
- The species name is extracted.
- Polynomial coefficients for high and low-temperature ranges are read and parsed from the subsequent lines.

4. Calculating Thermodynamic Properties:
- The Specific heat, enthalpy, and entropy functions are called to calculate Cp, H, and S for the given temperature range.
- These functions use the parsed polynomial coefficients and the specified temperature to compute the properties.

5. Molecular Weight Calculation:
- The CalMolWt function calculates the molecular weight of the species based on its chemical formula.

6. Plotting and Saving Results:
- The calculated thermodynamic properties are plotted against temperature.
- Separate figures are created for each property (Cp, H, S) and saved as images.

7. Output and File Handling:
- The molecular weights are written to a text file.
- A directory is created for each species, and the plots are saved in the respective directories.

This structured approach ensures that the code systematically reads, processes, calculates, and visualizes the thermodynamic data for various species, providing comprehensive insights into their properties over a specified temperature range.


Results

1. Molecular Weight of Species

2. Thermodynamic Property plots for O₂

Enthalpy for O2

Entropy for O2

Specific Heat For O2

3. Thermodynamic Property Plots for N₂

Enthalpy for N2

Entropy for N2

Specific Heat For N2

4. Thermodynamic Property Plots for CO₂

Enthalpy for CO2

Entropy for CO2

Specific Heat For CO2 Observations

5. Plots

Observations:

By analyzing the thermodynamic property plots of common gases like CO₂, O₂, and N₂, we can observe that properties such as enthalpy, entropy, and specific heat vary with temperature even under the ideal gas assumption.

For N2:

The diatomic nitrogen (N₂) shows a nearly linear variation in enthalpy with temperature. The entropy and specific heat, on the other hand, exhibit a somewhat logarithmic behavior. Overall, the thermodynamic properties of N2 continue to increase with rising temperature, with specific heat showing only slight variation above 2000 K. This steady increase is characteristic of diatomic gases as their vibrational modes become more accessible at higher temperatures.

For O2:

Diatomic oxygen (O₂) demonstrates a similar behavior to nitrogen in terms of enthalpy and entropy. However, the specific heat of O2 increases steeply up to 4000 K before starting to decrease. This behavior is due to the dissociation of O2 molecules into individual oxygen atoms (O) at higher temperatures. The energy required to break the O=O bond becomes available at these elevated temperatures, leading to a decrease in the specific heat capacity as the energy goes into breaking bonds rather than increasing kinetic energy.

For CO2:

Carbon dioxide (CO2) shows a trend in its thermodynamic properties similar to that of oxygen. CO2 also undergoes dissociation at temperatures above 4000 K, resulting in a decrease in specific heat. However, the variation in specific heat for CO2 is not as pronounced as it is for O2. This is likely due to the fact that CO₂, being a linear triatomic molecule, has more vibrational modes available at lower temperatures compared to diatomic molecules. Thus, its specific heat increases more gradually before the onset of dissociation.

These observations highlight the complexity and variability of thermodynamic properties with temperature, emphasizing the importance of accurate data and models in predicting the behavior of gases in high-temperature environments such as Combustion Research.

References:

1. https://ntrs.nasa.gov/citations/19950013764
2. https://www.osti.gov/servlets/purl/7073290
3. https://janaf.nist.gov/
4. https://publications.anl.gov/anlpubs/2005/07/53802.pdf
5. https://in.mathworks.com/matlabcentral/answers/196582-best-way-to-parse-text-file