NASA Thermodynamic File Parsing
Objective
1. To read and extract data present in thermo.dat data file.
2. To extract the 14 coeffecients of the given species.
3. To calculate enthalpy,entropy and specific heats using the mathematical relations.
4. To find the molecular weight of the given species.
5. To plot the specific heat,entropy and enthalpy relations for all species provided in the data.
6. To generate and save these properties and plots for each species in an appropriate directory.
Theory
Parsing is the process of analyzing text made of a sequence of tokens to determine its grammatical structure with respect to a given (more or less) formal grammar. The parser then builds a data structure based on the tokens. The process of "parsing" can also be described as reading in a data stream of some sort and building an in-memory model or representation of the semantic content of that data, in order to facilitate performing some kind of transformation on the data.
Here we are using NASA polynomial format for CHEMKIN II.The Chemkin file format is based very closely on NASA (1971) file format.
To take a look at the NASA Thermodynamic file click here NASA Thermodynamic data file
Code
Reading the thermodynamic data file
'fopen' command is used to open the file and the permission is set to read to read the file
'fgetl' command is used to jump from one line to another preffered line
ex. :- file_id = fopen(file,permission);
%reading the file and extracting the temperature values
f1 = fopen('thermo.dat','r');
fgetl(f1);Extracting the global temperature ranges
'strsplit' command is used to split the line containing the global temperature values based on the spaces (" ") .'str2num' command is used to convert the cell values of the temperature data extracted into numerical value so that it can be used for further calculations.
tdata = fgetl(f1);
t = strsplit(tdata,' ');
global_low_temp = str2num(t{2});
global_mid_temp = str2num(t{3});
global_high_temp = str2num(t{4});Calculating the number of lines
'while' loop and 'feof' command is used to read through the data file while seeking for the end. 'if' and 'break' command is used to sought the end of the file and match it for END ,once the condition is fulfilled the loop breaks. The total no. of lines are stored in the 'line' variable.
while ~feof(f1) %main loop to go for all substances
line = fgetl(f1);
a = sscanf(line,'%c');
%checking for end
c = a(1:3);
d = 'END'
if c == d;
break
end
line1 = line1 + 1;
end
Calculating the number of species
The total number of species is calculated using the following formula
%calculating the no. of species in the thermodynamics data file
no_species = (line1-5)/4;
Main body of the program
The main body of the program houses a loop which runs for all the 53 species whilst collecting all the relevant data and executing all the important functions. The functions used in the loop are :
1. 'get_coeff' is used to extract the local temperature ranges for the given species from the file, it extracts the name of each of the species and it also extracts all the 14 coefficients for each of the given specie.
2. A temperature array is created after extracting the local temperature ranges.
3. Temperature limits are also defined.
4. 'get_calc' function is used to calculate specific heat, entropy and enthalpy for all the species.
5. 'mol_wt' function is used to calculate the molecular weight of all the elements.
6. The molecular weight of each species is displayed.
7. 'plots' function is used to create a seperate folder for each specie and plot the various relations of cp,entropy ,enthalpy and temperature and dump it in their seprate folders.
for j = 1:no_species
[coeff,tempr,elem] = get_coeff(f1);
h_coeff = coeff(1:7);
l_coeff = coeff(8:14);
t1 = tempr(1);
t2 = tempr(2);
t3 = tempr(3);
temp = linspace(t1,t2,300);
m_temp = tempr(3);
[cp,h,s] = get_calc(h_coeff,l_coeff,temp,m_temp);
Cp = cp;
H = h;
S = s;
element = elem;
molecular_weight(j) = mol_wt(element);
fprintf('The Molecular weight of %s is %4.3f', element,molecular_weight(j))
molar_mass = molecular_weight(j);
plots(Cp,H,S,temp,element,molar_mass);
end
Functions
Function to extract various parameters : get_coeff function
This function extracts the name for each of the species , extracts their local temperature ranges and it also the extracts the 14 coefficients all the elements.
1. species variable stores the name of each specie.
2. 'If' loop is used to extract the local temperature ranges using the 'strfind' command to find the location of the G in the data file.
3. 'fgetl' command is used to change the lines.
4. 'strfind' command is used to find the location of E in the coefficients which is then use to point out the entire value and extract it.
5. 'str2double' command is used to convert the data into integer form so that it can be used for calculations.
function[coeff,tempr,elem] = get_coeff(f1)
name_line = fgetl(f1);
species = name_line(1:6);
templ = strfind(name_line,'G');
if templ(1)<45
temp1 = name_line(templ(1)+30:templ(1)+36);
temp1 = str2double(temp1);
temp2 = name_line(templ(1)+39:templ(1)+46);
temp2 = str2double(temp2);
temp3 = name_line(templ(1)+49:templ(1)+56);
temp3 = str2double(temp3);
else
temp1 = name_line(templ+4:templ+10);
temp1 = str2double(temp1);
temp2 = name_line(templ+13:templ+20);
temp2 = str2double(temp2);
temp3 = name_line(templ+23:templ+30);
temp3 = str2double(temp3);
end
A = fgetl(f1);
B = fgetl(f1);
C = fgetl(f1);
%finding location of the E in the above strings
a = strfind(A,'E');
b = strfind(B,'E');
c = strfind(C,'E');
%getting all the 14 coefficients
a1 = A(1:a(1)+3);
a1 = str2double(a1);
a2 = A(a(1)+4:a(2)+3);
a2 = str2double(a2);
a3 = A(a(2)+4:a(3)+3);
a3 = str2double(a3);
a4 = A(a(3)+4:a(4)+3);
a4 = str2double(a4);
a5 = A(a(4)+4:a(5)+3);
a5 = str2double(a5);
a6 = B(1:b(1)+3);
a6 = str2double(a6);
a7 = B(b(1)+4:b(2)+3);
a7 = str2double(a7);
a8 = B(b(2)+4:b(3)+3);
a8 = str2double(a8);
a9 = B(b(3)+4:b(4)+3);
a9 = str2double(a9);
a10 = B(b(4)+4:b(5)+3);
a10 = str2double(a10);
a11 = C(1:c(1)+3);
a11 = str2double(a11);
a12 = C(c(1)+4:c(2)+3);
a12 = str2double(a12);
a13 = C(c(2)+4:c(3)+3);
a13 = str2double(a13);
a14 = C(c(3)+4:c(4)+3);
a14 = str2double(a14);
coeff = [a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14];
tempr = [temp1 temp2 temp3];
elem = species;
end
Function to calculate molecular weight : mol_wt function
This function takes the specie name as the input and calculates their weight.This function is explained as follows.
1. An array 'p' is created to store the name of single elements then another array 'at_wt' is created to store the weights of the single elements which were used earlier .The weights are stored at the same index values of the elemnts decribed earlier.
2. Simple for loop is run for times equal to length of 'specie' array which is provided at the input. Another nested for loop is run for times equal to length of 'p'.
3. In nested loop each letter of the given spiecies name is compared with element name stored in 'p' then using if statement it is found out if its true or false. If its true molecular weight gets added and nested loop breaks ,if it is found false it repeats the loop. 'str2double' command is used to find out if spiecies contains any numeric value. 'if' statement is used again to check if value of n is greater than 1. If its true then molecular weight get adjusted by following formula mentioned in code. If false for loop repeats itself again.
function [w] = mol_wt(specie)
p = ['H','C','N','O','A'];
at_wt = [1 ,12 ,14 ,16, 40];
w=0;
%calculating molecular weight of the species
for i = 1:length(specie)
for j = 1:length(p)
if(strcmp(specie(i),p(j))==1)
w=w+at_wt(j);
J=j;
break;
end
end
q=str2double(specie(i));
if (q>1)
w=w+(at_wt(J).*(q))-at_wt(J);
end
end
endFunction to calculate the specific heat, entropy and enthalpy : get_ calc function
This funtion is used to calculate the value of specific heat, entropy and enthalpy for each specie. The function is explained below:
1. value of Universal gas constant is defined at the begining which is later to be used in the calculations
2. The following formulas are used to find the specific heat,enthalpy and entropy of the given specie.
3. A for loop which runs from intial value of 1 to the length of entire temperature range array is used.An 'if else' loops is used to diffrentiate the coefficients as high temperature and low temperature coefficients based on the temperature limit defined. Various values of the extracted coefficients using the 'get_coeff' is given to the function as an input.
function [cp,h,s] = get_calc(ht,lt,t,m_temp)
R = 8.314; %unit - J/mol*k
for k = 1:length(t)
if t(k)>m_temp
cp(k) = R*(ht(1) + (ht(2)*(t(k)) + (ht(3)*(t(k)^2)) + (ht(4)*((t(k))^3)) + (ht(5)*((t(k))^4))));
h(k) = R*((ht(1)*(t(k))) + (ht(2)*((t(k))^2)/2) + (ht(3)*((t(k)^3)/3) + (ht(4)*((t(k))^4)/4) + (ht(5)*((t(k))^5)/5) + ht(6)));
s(k) = R*((ht(1)*log(t(k))) + (ht(2)*(t(k))) + (ht(3)*((t(k))^2)/2) + (ht(4)*((t(k))^3)/3) + (ht(5)*((t(k))^4/4) + ht(7)));
else
cp(k) = R*(lt(1) + (lt(2)*(t(k)) + (lt(3)*(t(k)^2)) + (lt(4)*((t(k))^3)) + (lt(5)*((t(k))^4))));
h(k) = R*((lt(1)*(t(k))) + (lt(2)*((t(k))^2)/2) + (lt(3)*((t(k)^3)/3) + (lt(4)*((t(k))^4)/4) + (ht(5)*((t(k))^5)/5) + ht(6)));
s(k) = R*((lt(1)*log(t(k))) + (lt(2)*(t(k))) + (lt(3)*((t(k))^2)/2) + (lt(4)*((t(k))^3)/3) + (lt(5)*((t(k))^4/4) + lt(7)));
end
end
end
Function to generate the plots for various relations: get_ plot function
This function is used to plot various graphs between various relations and store them in a seperate file.
1. This function takes entropy,specific heat , enthalpy and temperature as inputs to generate the plots for the given species.
2. Two strings "properties of' and element are merged using the 'strcat' command and the mkdir command is used to create the directory with that name.
3. The graphs are plotted for various relations using the 'plot' command and 'saveas' command is used to save the plots whilst giving the designated name as per the relations it's ploting the graph for.'Fullfile' cmmand is used to set the location of the file. After saving the files the directory is changed to the orginal directory using the 'cd' command.
4. The same process is repeated for each relation and all the 3 relational graphs are obtained.
5. Molecular weights of the given species is stored in a seperate text file using the 'fopen' command and later writing the value using the write permission.
function plots(cp,h,s,temp,element,molar_mass)
data=strcat('properties of',element);
mkdir(data);
file = pwd;
%plotting and saving data
figure(1)
plot(temp,cp);
data1=strcat(element,'Specific Heat vs temperature');
legend(data1,'location','northoutside');
xlabel('Temperature(K)');
ylabel('Specific Heat(cp) [J/mol*k]');
f=fullfile(data);
cd(f);
saveas(figure(1),[pwd '/specific_heat_vs_temp_plot.jpg']);
cd(file);
figure(2)
plot(temp,h);
data2 = strcat(element,'Enthalpy vs Temperature');
legend(data2,'location','northoutside');
xlabel('Temperature(K)');
ylabel('Enthalpy(h) [J/mol]');
cd(f);
saveas(figure(2),[pwd '/Enthalpy_vs_Temp_plot.jpg']);
cd(file);
figure(3)
plot(temp,s);
data3 = strcat(element,'Entropy vs Temperature');
legend(data3,'location','northoutside');
xlabel('temperature(K)');
ylabel('Entropy(s) [J/mol*k]');
cd(f)
saveas(figure(3),[pwd '/Entropy_vs_temp_plot.jpg']);
cd(file)
%saving the molecular weight of the species in the file
cd(f)
wt1 = fopen('Molecular Mass.txt','wt');
fprintf(wt1,'%s %s %f \n','Molecular Weight of',element,molar_mass);
fclose(wt1);
cd(file)
end
Complete Code synopsis
%NASA thermodynamic file parsing
close all
clear all
clc
%reading the file and extracting the temperature values
f1 = fopen('thermo.dat','r');
fgetl(f1);
tdata = fgetl(f1);
t = strsplit(tdata,' ');
global_low_temp = str2num(t{2});
global_mid_temp = str2num(t{3});
global_high_temp = str2num(t{4});
%skipping the next few comment line
for i = 1:3
fgetl(f1);
end
frewind(f1);
line1 = 0;
while ~feof(f1) %main loop to go for all substances
line = fgetl(f1);
a = sscanf(line,'%c');
%checking for end
c = a(1:3);
d = 'END'
if c == d;
break
end
line1 = line1 + 1;
end
%calculating the no. of species in the thermodynamics data file
no_species = (line1-5)/4;
frewind(f1);
for i = 1:5 %skiping unwanted lines
fgetl(f1);
end
for j = 1:no_species
[coeff,tempr,elem] = get_coeff(f1);
h_coeff = coeff(1:7);
l_coeff = coeff(8:14);
t1 = tempr(1);
t2 = tempr(2);
t3 = tempr(3);
temp = linspace(t1,t2,300);
m_temp = tempr(3);
[cp,h,s] = get_calc(h_coeff,l_coeff,temp,m_temp);
Cp = cp;
H = h;
S = s;
element = elem;
molecular_weight(j) = mol_wt(element);
fprintf('The Molecular weight of %s is %4.3f', element,molecular_weight(j))
molar_mass = molecular_weight(j);
plots(Cp,H,S,temp,element,molar_mass);
end
function[coeff,tempr,elem] = get_coeff(f1)
name_line = fgetl(f1);
species = name_line(1:6);
templ = strfind(name_line,'G');
if templ(1)<45
temp1 = name_line(templ(1)+30:templ(1)+36);
temp1 = str2double(temp1);
temp2 = name_line(templ(1)+39:templ(1)+46);
temp2 = str2double(temp2);
temp3 = name_line(templ(1)+49:templ(1)+56);
temp3 = str2double(temp3);
else
temp1 = name_line(templ+4:templ+10);
temp1 = str2double(temp1);
temp2 = name_line(templ+13:templ+20);
temp2 = str2double(temp2);
temp3 = name_line(templ+23:templ+30);
temp3 = str2double(temp3);
end
A = fgetl(f1);
B = fgetl(f1);
C = fgetl(f1);
%finding location of the E in the above strings
a = strfind(A,'E');
b = strfind(B,'E');
c = strfind(C,'E');
%getting all the 14 coefficients
a1 = A(1:a(1)+3);
a1 = str2double(a1);
a2 = A(a(1)+4:a(2)+3);
a2 = str2double(a2);
a3 = A(a(2)+4:a(3)+3);
a3 = str2double(a3);
a4 = A(a(3)+4:a(4)+3);
a4 = str2double(a4);
a5 = A(a(4)+4:a(5)+3);
a5 = str2double(a5);
a6 = B(1:b(1)+3);
a6 = str2double(a6);
a7 = B(b(1)+4:b(2)+3);
a7 = str2double(a7);
a8 = B(b(2)+4:b(3)+3);
a8 = str2double(a8);
a9 = B(b(3)+4:b(4)+3);
a9 = str2double(a9);
a10 = B(b(4)+4:b(5)+3);
a10 = str2double(a10);
a11 = C(1:c(1)+3);
a11 = str2double(a11);
a12 = C(c(1)+4:c(2)+3);
a12 = str2double(a12);
a13 = C(c(2)+4:c(3)+3);
a13 = str2double(a13);
a14 = C(c(3)+4:c(4)+3);
a14 = str2double(a14);
coeff = [a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14];
tempr = [temp1 temp2 temp3];
elem = species;
end
function [w] = mol_wt(specie)
p = ['H','C','N','O','A'];
at_wt = [1 ,12 ,14 ,16, 40];
w=0;
%calculating molecular weight of the species
for i = 1:length(specie)
for j = 1:length(p)
if(strcmp(specie(i),p(j))==1)
w=w+at_wt(j);
J=j;
break;
end
end
q=str2double(specie(i));
if (q>1)
w=w+(at_wt(J).*(q))-at_wt(J);
end
end
endfunction plots(cp,h,s,temp,element,molar_mass)
data=strcat('properties of',element);
mkdir(data);
file = pwd;
%plotting and saving data
figure(1)
plot(temp,cp);
data1=strcat(element,'Specific Heat vs temperature');
legend(data1,'location','northoutside');
xlabel('Temperature(K)');
ylabel('Specific Heat(cp) [J/mol*k]');
f=fullfile(data);
cd(f);
saveas(figure(1),[pwd '/specific_heat_vs_temp_plot.jpg']);
cd(file);
figure(2)
plot(temp,h);
data2 = strcat(element,'Enthalpy vs Temperature');
legend(data2,'location','northoutside');
xlabel('Temperature(K)');
ylabel('Enthalpy(h) [J/mol]');
cd(f);
saveas(figure(2),[pwd '/Enthalpy_vs_Temp_plot.jpg']);
cd(file);
figure(3)
plot(temp,s);
data3 = strcat(element,'Entropy vs Temperature');
legend(data3,'location','northoutside');
xlabel('temperature(K)');
ylabel('Entropy(s) [J/mol*k]');
cd(f)
saveas(figure(3),[pwd '/Entropy_vs_temp_plot.jpg']);
cd(file)
%saving the molecular weight of the species in the file
cd(f)
wt1 = fopen('Molecular Mass.txt','wt');
fprintf(wt1,'%s %s %f \n','Molecular Weight of',element,molar_mass);
fclose(wt1);
cd(file)
end
function [cp,h,s] = get_calc(ht,lt,t,m_temp)
R = 8.314; %unit - J/mol*k
for k = 1:length(t)
if t(k)>m_temp
cp(k) = R*(ht(1) + (ht(2)*(t(k)) + (ht(3)*(t(k)^2)) + (ht(4)*((t(k))^3)) + (ht(5)*((t(k))^4))));
h(k) = R*((ht(1)*(t(k))) + (ht(2)*((t(k))^2)/2) + (ht(3)*((t(k)^3)/3) + (ht(4)*((t(k))^4)/4) + (ht(5)*((t(k))^5)/5) + ht(6)));
s(k) = R*((ht(1)*log(t(k))) + (ht(2)*(t(k))) + (ht(3)*((t(k))^2)/2) + (ht(4)*((t(k))^3)/3) + (ht(5)*((t(k))^4/4) + ht(7)));
else
cp(k) = R*(lt(1) + (lt(2)*(t(k)) + (lt(3)*(t(k)^2)) + (lt(4)*((t(k))^3)) + (lt(5)*((t(k))^4))));
h(k) = R*((lt(1)*(t(k))) + (lt(2)*((t(k))^2)/2) + (lt(3)*((t(k)^3)/3) + (lt(4)*((t(k))^4)/4) + (ht(5)*((t(k))^5)/5) + ht(6)));
s(k) = R*((lt(1)*log(t(k))) + (lt(2)*(t(k))) + (lt(3)*((t(k))^2)/2) + (lt(4)*((t(k))^3)/3) + (lt(5)*((t(k))^4/4) + lt(7)));
end
end
end
Plots Generated
Plots for O2
Plots for N2
Plots for CO2
Screenshot of the folders in which the plots and weights are stored
Errors encountered
While running the program initially the program was unable to plot the graphs for 2 elemnts viz. HCN & CH2OH ,with further inspection and analysis it was found that these 2 elements had two values of G in the line in which their local temperature has been defined.Therefore it was facing issues to extract the temperature ranges and hence it was unable to plot the graphs.
The screen shot of the error can be seen below:
Error Rectification
The above error was eliminated by using a 'if else' condition statement inside the 'get_coeff' function where it extracts the temperatures.The 'if' command looks out for values lessser than the conventional locational value of G i.e 45 and as soon as it encounters the case which fullfils the conditions a sperate set of commands is excuted which extracts the temperature ranges for CH2OH and HCN, if no such condition is encountered it executes the 'else' conditional statement and therefore program is able to extract each and every temperature ranges for the given species.