Parsing NASA's thermodynamic data file using MATLAB

Aim:

To parse and create value out of NASA's thermodynamic data

Objectives:

1. Extract the 14 coefficients and calculate specific heat (Cp), enthalpy (H), and entropy (S) for all the species in the given data
2. Calculate the molecular weight of each species and display it in the command window
3. Plot the Cp, H, and S with respect to the local temperature range unique to each species
4. Save the plots as images, with each species having its own folder

Theory:

File parsing is the process of building an in-memory model or representation of the data in the respective file, in order to analyse and use the data to perform some kind of transformation. Parsing can be used to both read and write data from/to a file.

NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S. They have also documented the co-efficients that are required to evaluate these polynomials for 53 species, as can be referred to in the file above. This project is about creating sense out of all of that information, as listed in the objectives.

Governing Equations:

`C_p = R(a_1 + a_2T + a_3T^2 + a_4T^3 + a_5T^4)`

`H = RT(a_1 + (a_2T)/2 + (a_3T^2)/3+ (a_4T^3)/4+ (a_5T^4)/5+a_6/T)`

`S = R(a_1lnT + a_2T+ (a_3T^2)/2 + (a_4T^3)/3+ (a_5T^4)/4 +a_7)`

Code for Objective A:

clear all; close all; clc

f2 = fopen('THERMO.dat','r');
%------------------Parsing the header info------------------
for i = 1:2
fgetl(f2);
end

%extracting global temperatures
global_temperatures = str2num(ans);
global_low_temp = global_temperatures(1);
global_mid_temp = global_temperatures(2);
global_high_temp = global_temperatures(3);

%skipping comment section (lines 3 to 5)
for i = 1:3
    fgetl(f2);
end

%---------------Parsing the species block-------------------
for i = 1:53
species_info = fgetl(f2); %to obtain the entire string of information
species_info = strsplit(species_info,' '); %to split the string into separate elements via a cell array
species_chemical_formula = species_info{1}; %first element of the string

%extracting the local temp ranges that is unique to each species
local_low_temp = str2double(species_info{end - 3});
local_high_temp = str2double(species_info{end - 2});
local_mid_temp = str2double(species_info{end - 1});

%----extracting coefficients a1 to a5-----
coeff_line_1 = fgetl(f2); 
a = strfind(coeff_line_1,'E');
a1 = str2double(coeff_line_1(1:a(1)+3));
a2 = str2double(coeff_line_1(a(1)+4:a(2)+3));
a3 = str2double(coeff_line_1(a(2)+4:a(3)+3));
a4 = str2double(coeff_line_1(a(3)+4:a(4)+3));
a5 = str2double(coeff_line_1(a(4)+4:a(5)+3));

%----extracting coefficients a6 to a10-----
coeff_line_2 = fgetl(f2); 
a = strfind(coeff_line_2,'E');
a6 = str2double(coeff_line_2(1:a(1)+3));
a7 = str2double(coeff_line_2(a(1)+4:a(2)+3));
a8 = str2double(coeff_line_2(a(2)+4:a(3)+3));
a9 = str2double(coeff_line_2(a(3)+4:a(4)+3));
a10 = str2double(coeff_line_2(a(4)+4:a(5)+3));

%----extracting coefficients a11 to a14-----
coeff_line_3 = fgetl(f2);
a = strfind(coeff_line_3,'E');
a11 = str2double(coeff_line_3(1:a(1)+3));
a12 = str2double(coeff_line_3(a(1)+4:a(2)+3));
a13 = str2double(coeff_line_3(a(2)+4:a(3)+3));
a14 = str2double(coeff_line_3(a(3)+4:a(4)+3));

T = linspace(local_low_temp,local_high_temp,1000); %defining temp range
R = 8.3145; %universal gas constant

for j = 1:length(T)
if T(j) > local_mid_temp
    Cp(j) = R*(a1 + a2.*T(j) + a3.*T(j).^2 + a4.*T(j).^3 + a5.*T(j).^4); %high temp coeffs
    H(j) = R.*T(j).*(a1 + (a2.*T(j))./2 + (a3.*T(j).^2)./3 + (a4.*T(j).^3)./4 + (a5.*T(j).^4)./5 + a6./T(j));
    S(j) = R*(a1.*log(T(j)) + a2.*T(j) + (a3.*T(j).^2)./2 + (a4.*T(j).^3)./3 + (a5.*T(j).^4)./4 + a7);
else 
    Cp(j) = R*(a8 + a9.*T(j) + a10.*T(j).^2 + a11.*T(j).^3 + a12.*T(j).^4); %low temp coeffs
    H(j) = R.*T(j).*(a8 + (a9.*T(j))./2 + (a10.*T(j).^2)./3 + (a11.*T(j).^3)./4 + (a12.*T(j).^4)./5 + a13./T(j));
    S(j) = R*(a8.*log(T(j)) + a9.*T(j) + (a10.*T(j).^2)./2 + (a11.*T(j).^3)./3 + (a12.*T(j).^4)./4 + a14);
end
end

The code under 'Parsing the species block' is what is relevant to meet Objective A, and will be the discussed part.

It is key to realize that the same repetitive procedure will be used to parse through the data for all 53 species, so a code that works for all must be written, and hence the use of the for-loop. Following that, all the data can only be extracted as a string at first, using commands such as strsplit and strfind, which then needs to be converted into text or numbers accordingly.

Referring to the data file, the 14 coefficients given for each species are distributed in 3 lines, with 5 coefficients each in the first and second line. The first 7 coefficients are for the high-temperature range and the next 7 coefficients are for the low-temperature range. Each species had unique local temperatures, which also needed to be extracted and stored in order to create an array of the temperature range.

Another for-loop needed to be incorporated to scan through each of the temperature values in the range that was determined. An if-else statement was also necessary to use the appropriate coefficients for the low and high temperature values in order to calculate the accurate values of Cp, H, and S.

Code for Objective B:

%----------Calculating molecular weight of each species----------

atoms = ['C','H','O','N','Ar'];
atomic_mass = [12,1,16,14,40];
M = 0;
for m = 1:length(species_chemical_formula) %to scan through each constituent of the species' chemical formula
    for n = 1:length(atoms) %to scan through all the possible atoms defined above
          if(strcmp(species_chemical_formula(m),atoms(n))) %comparing whether each constituent matches with the atom
              M = M + atomic_mass(n); %calculating molecular weight via summation
              store = n; %storing the value of the position of n in the atoms array; this is very important
          else
          end
    end
    number = str2num(species_chemical_formula(m)); 
    if (number > 1)
        M = M + atomic_mass(store)*(number-1); %calculating molecular weight for the species that have more than one atom of an element
    else
    end
end

%-----displaying molecular weight in command window-----

fprintf('Molecular weight of %s: ',species_chemical_formula);
fprintf('%f',M);
disp(' ')

Calculating the molecular weight of each species meant that the first step was to identify the constituents that were present in all of the 53 species. The atoms array showcases those constituents. Once again, it was to be identified that the molecular weight of each species will be calculated in the same manner, and hence the use of the for-loop.

What was crucial to realize though was that, every constituent was to be scanned through for each species. Since it was already determined that only 5 types of constituents existed in the entire data, the scanned constituents needed to be compared to the constituents defined in the atoms array using the strcmp command to see if there was a match. This is why a nested-for loop was needed. If there was a match, the molecular weight was to be calculated, and the process continued for all the atoms in a single species before moving on to the next species. This is why an if-else statement was needed. 

So far, the text-only data for each species such as CH or OH were accounted for, but there are multiple atoms of a specific species for several species, such as C3H7 and H2O2. Another if-else statement was needed to account for those multiple atoms.

Code for Objective C and D:

%-----creating a separate directory to store plots and molecular weight of all species------
cd('C:UsersParth MaheshwariDesktopSKILL LYNC MC-HEVMATLAB FOR MECHANICAL ENGINEERS (1)tempFile Parsing Project');
mkdir(species_chemical_formula);
cd(species_chemical_formula);

%------writing molecular weight to a text file-----
MW_file = fopen('Molecular Weight.txt','w');
fprintf(MW_file,'Molecular Weight of %s is %d',species_chemical_formula,M);
fclose(MW_file);

%-------------Plotting graphs for each species-------------
%Temp vs Cp
figure;
subplot(3,1,1); %syntax represents 3 rows, 1 column, 1st position
plot(T,Cp,'linew',1.5,'color','r');
xlabel('Temperature (K)');
ylabel('Specific Heat (KJ)');
title(sprintf('Specific Heat vs Temperature of %s',species_chemical_formula));

%Temp vs Enthalpy
subplot(3,1,2);
plot(T,H,'linew',1.5,'color','b');
xlabel('Temperature (K)');
ylabel('Enthalpy (KJ/mol)');
title(sprintf('Enthalpy vs Temperature of %s',species_chemical_formula));

%Temp vs Entropy
subplot(3,1,3);
plot(T,S,'linew',1.5,'color','k');
xlabel('Temperature (K)');
ylabel('Entropy (KJ/mol-K)');
title(sprintf('Entropy vs Temperature of %s',species_chemical_formula));

saveas(figure(i),'Cp,H,S vs Temp.jpg')

%changing directory to original
cd('C:UsersParth MaheshwariDesktopSKILL LYNC MC-HEVMATLAB FOR MECHANICAL ENGINEERS (1)temp')

Before plotting Cp, H, and S with respect to temperature, the directory was changed to a new folder where all the output was to be saved. This was done using the cd command. Within the new directory, the mkdir command was used to create a new folder for each species' output plot and molecular weight. The cd command was used once again here to change the directory to the specific species for output saving purposes. These steps ensured that when the plots were generated and saved, each species' plot would be in its own folder.

                                                            Output: All folders saved in a separate directory named "File Parsing Project"

                                                                                                                      Output plots of `O_2, N_2, CO_2`:

                                                                                      Output: Molecular weight of each species displayed in command window