Project 1 - Parsing NASA thermodynamic data
Objective
1) To Read ,Extract and operate Data Present in thermo30.data
2) To extract the 14 co-efficients for a given species.
3)To calculate the enthalpy, entropy and specific heats using the NASA polynomials.
4)To calculate the molecular weight of each species.
5) To Plot the Cp, Enthalpy and Entropy of all species present in Thermo Data. 6) To save plots of particular species in a appropriate directory.
7) To Generate and save table containing molecular weight and Species Name.
Theory
Parsing is the process of analysing a string of symbols, either in natural language, computer languages or data structures, conforming to the rules of a formal grammar.Within computational linguistics the term is used to refer to the formal analysis by a computer of a sentence or other string of words into its constituents, resulting in a parse tree showing their syntactic relation to each other, which may also contain semantic and other information.
[1] In the following Program we are using NASA Polynomial format for CHEMKIN-II.The Chemkin thermo file format is based very closely on the NASA (1971) file format [2][3].To read more about chemkin thermo file format click Read about Chemkin Format.
The view NASA Thermodynamic Data File[4] Click NASA Thermodynamic Data
Code synthesis
Reading thermodynamic file
The file can be opened using fopen command [5] and permission is set to read
filelD = fopen(filename,permission)
e.g f
l=fopen('thermo30.dat.txt7r.)
Calculating No of lines
Using for loop running infinite times, the program reads file line by line by using fgetl command[6]. Then it compares the line with term "END' to find out the end of file using strcmp command[7].If it returns true then loop breaks and value of no of lines is equal to count of loop at break[8].
e.g
for i=1:Inf
line= fgetl(f1); flag=strcmp
(line,'END'); if(flag==1)
no_of_lines=i; break;
end
end
Calculating No of species
present It can be calculated by simple formula
no_of_species=(no_of_lines-6)/4;
Finding global temperature range
Finding Global Temperature Range To find global temperature range, first file pointer is set to initial line by using frewind command [9].Then first line is dumped and second line is read and stored in String. Then this string is broken from locations of space(' ')by using strsplit command[10] and then correct strings from broken strings are selected and converted it to numerical value by str2double command[11].
e.g
%finding Global Temperature Range
frewind(fl); % to bring file pointer to intial position of file
fgetl(f1); %reading the first line
data = fgetl(fl); %reading the second
line split = strsplit(data,"); % separating string by space
for i = 1:length(split)
gt(i) = str2double(split{i}); %convrsion of string charaters to numbers
end
global_min = gt(2); %global min temperature
limit = gt(3); %global med temperature
global_high = gt(4); %global max temperature
Creating a main body of program
Main body of program runs in for loop ,for times equal to no of species present in data.This loop contains various functions to get diffferent data which are as follows
a) get_coff is function to evaluate all 14 coefficients of a particular species.Then first 7 coefficient gets stored in h_coff and last 7 coefficient get stored in 1_coff.
b) A time array is created using linspace command.
c)get_cp, get_h, get_s are the functions to evaluate Cp , enthalpy and entropy respectively for a particular spieces.
d)get_element is a function to get spieces Name for a particular spieces.
e)get_mw is a function to get spieces molecular weight for a particular spieces.
f)get_plots is a function to plot spieces data and save result figures in jpg format for a particular spieces.
g)Chemical_formula is an array to store all spieces name. Transpose of Chemical_formula and molecular_wt is used to create a table by using table command[12] and table in written in xlsx format by using writetable command [13].
for k=i:no_of_species
str=k;
coff=get_coff(fl,str); %function to get 14 coefficients
h_coff=coff(1:7); % high Temperature coefficients
l_coff=coff(8:14); % low Temperature coefficient
temp=linspace(300,3500,320); % this temperature range is a best fit between local and global temperature
cp=get_cp(h_coff,l_coff,temp,limit); % function to attain cp value for particular Spieces
H=get_h(h_coff,l_coff,temp,limit); % function to attain Enthapy value for particular Spieces
S=get_s(h_coff,l_coff,temp,limit); % function to attain Entropy value for particular Spieces
element=get_element(fl,str); % function to attain spiecies Name value for particular Spieces molecular_wt(k)=get_mw(element) ; % function to attain spiecies Molecular Weight value for particular Spieces
get_plots(cp,H,S,temp,element); %function to attain Plots and saving plots to relevant location for particular Speices
chemical_formula(k)={get_element(fl,str)}; %storing spieces Name in a cell array
end
%generating table to store results of molecular weight
Spieces_Molecular_weight = table(chemical_formula', molecular_wt'); %generating table
Spieces_Molecular_weight.Properties.VariableNames = {'Spieces' Molecular_weight' }; %Naming table variables
writetable(Spieces_Molecular_weight,'Spieces_Molecular_weight.xlsx'); % saving table in xlsx format.
Creating get_coeff function
This function can be understood by following statements
a) 'i' indicates line number where coefficient data of particular species is stored.
b) File pointer is set to initial line by using frewind command.
c) Using textscan command[14] data can be read and stored from ith line as shown below. Then using fgetl command data of next two lines is stored.
d) Since textscan creates a cell arrayof lxl , the cell array 'd' is converted into char string using char command and since is cell array '{}' is used.
e)Location of exponent E is found out using strfind command[15],then first coefficient is stored in string format and then it is converted to number by using str2double command. Same procedure is caried out to find out all remaining coefficient.Then an array is created to store all 14 coefficients.
function[coffl=get_coff(f1,str)
i=4*(str-1)+7; % getting spiecies Number
frewind(f1); % to bring file pointer to intial position of file
d = textscan(f1,'%s',1,'delimiter','n', 'headerlines',i-1); % to read line from i position
D=char(d{1}); %converting cell array to string
E=fgetl(f1); % Reading next line
F=fgetl(f1); % Reading next line
%Finding Location of Exponent in D,E,F
string a=strfind(D,'E');
b=strfind(E,'E');
c=strfind(F,'E');
%getting 14 coefficients
a1=D(1:a(1)+3);
a1=str2double(a1);
a2=D(a(1)+4:a(2)+3);
a2=str2double(a2);
a3=D(a(2)+4:a(3)+3);
a3=str2double(a3);
a4=D(a(3)+4:a(4)+3);
a4=str2double(a4);
a5=D(a(4)+4:a(5)+3);
a5=str2double(a5);
a6=E(1:b(1)+3);
a6=str2double(a6);
a7=E(b(1)+4:b(2)+3);
a7=str2double(a7);
a8=E(b(2)+4:b(3)+3);
a8=str2double(a8);
a9=E(b(3)+4:b(4)+3);
a9=str2double(a9);
a10=E(b(4)+4:b(5)+3);
a10=str2double(a10);
a11=F(1:c(1)+3);
a11=str2double(a11);
a12=F(c(1)+4:c(2)+3);
a12=str2double(a12);
a13=F(c(21+4:c(3)+31:
a13=str2double(a13);
al4=F(c(3)+4:c(4)+3);
al4=str2double(a14);
coff=[al a2 a3 a4 a5 a6 a7 a8 a9 al0 all a12 a13 a14];
Creating get_cp function
This function takes inputs as high temperature coefficients , low temperature coefficients, temprature array and limit to decide whether to use high cofficient or low coefficients. This functions contains a simple for loop to carry out operations at each temperature in temperature array. Cp can be calculated by using
Cp/R = al + a2 T + a3 TA2 + a4 TA3 + a5 TA4
then if condition is used to decide, whether to use high coefficient or low coefficients and it is decided by value stored in limit variable.
function [Cp]=get_cp(high_temp,low_temp,temp,limit)
R=8.314; %unit j/mol*K
% Calculating Cp for different temperatures
for i=1:length(temp)
terml = high_temp(1);
term2 = high_temp(2).*temp(i);
term3 = high_temp(3).xtemp(i)^2;
term4 = high_temp(4).*temp(i)^3;
term5 = high_temp(5).xtemp(i)^4;
term6 = low_temp(1);
term7 = low_temp(2).*temp(i);
term8 = low_temp(3).*temp(i)^2;
term9 = low_temp(4).*temp(i)^3;
term10 = low_temp(5).*temp(i)^4;
if temp(i)>=1imit
Cp(i)= R.*(term1+ term2 +term3+ term4+term5);
else
Cp(i)=R.*( term6 +term7 +term8 +term9 +term10);
end
end
end
Creating get_h function
This function takes inputs as high temperature coefficients , low temperature coefficients, temprature array and limit to decide whether to use high cofficient or low coefficients. This functions contains a simple for loop to carry out operations at each temperature in temperature array. Enthalpy can be calculated by using
H/RT = al + a2 T /2 + a3 T^2 /3 + a4 T^3 /4 + a5 T^4 /5 + a6/T
then if condition is used to decide, whether to use high coefficient or low coefficients and it is decided by value stored in limit variable.
function [H]=get_h(high_temp,low_temp,temp,limit)
R=8.314; %unit 7/mol*K
% Calculating enthalpy for different temperatures
for i=1:length(temp)
term1=high_temp(1).*temp(i);
term2=(high_temp(2).*temp(i)^2)/2;
term3=(high_temp(3).*(temp(i)^3))/3;
term4=(high_temp(4).*(temp(i)^4))/4;
term5=(high_temp(5).*(temp(i)^5))/5;
term6=high_temp(6);
term7=low_temp((1).*temp(i);
term8=(low_temp((2).*temp(i)^2)/2;
term9=(low_temp((3).*(temp(i)^3))/3;
term10=(low_temp((4).*(temp(i)^4))/4;
ter11=(low_temp((5).*(temp(i)^5))/5;
term12=low_temp(6);
if temp(i)>=limit
H(1)=R.*(term1+term2+term3+term4+term5+term6);
else
H(1)=R.*(term7+term8+term9+term10+term11+term12);
end
end
end
Creating get_s function
This function takes inputs as high temperature coefficients , low temperature coefficients, temprature array and limit to decide whether to use high cofficient or low coefficients. This functions contains a simple for loop to carry out operations at each temperature in temperature array. Entropy can be calculated by using
S/R = al InT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7
then if condition is used to decide, whether to use high coefficient or low coefficients and it is decided by value stored in limit variable.
function [S]=get_s(high_temp,low_temp,temp,limit)
R=8.314; %unit 7/mol*K
% Calculating Entropy for different temperatures
for i=1:length(temp)
term1=high_temp(1).*(log(temp(i)));
term2=high_temp(2).*temp(i);
term3=(high_temp(3).*(temp(i)^2))/2;
term4=(high_temp(4).*(temp(i)^3))/3;
term5=(high_temp(5).*(temp(i)^4))/4;
term6=high_temp(7);
term7=low_temp(1).*(log(temp(i)));
term8=low_temp(2).*temp(i);
term9=(low_temp(3).*(temp(i)^2))/2;
term10=(low_temp(4).*(temp(i)^3))/3;
term11=(low_temp(5).*(temp(i)^4))/4;
term12=high_temp(7);
if temp(i)>=1imit
S(i)=R.*(term1+term2+term3+term4+term5+term6);
else
S(1)=R.*(term7+term8+term9+term10+term11+term12);
end
end
end
Creating get_element function
This function can be understood by following statements:-
a) 'i' indicates line number where species name of particular species is stored.
b) File pointer is set to initial line of file by using frewind command.
c) Using textscan command data can be read and stored from ith line as shown below.
d) Since textscan creates a cell arrayof lxl , the cell array 'r' is converted into char string using char command and since is cell array'{}' is used.
e)Location of space C ')is found out using strfind command ,then string is stored in 'ele' by using simple transformation shown below in code.
function [ ele]=get_element (fl, str)
i=4*(str-1)+6; % getting spiecies Number
frewind(f1); % Setting file pointer to intial position of file
r= textscan(f1,'%s',1,'delimiter','n','headerlines',i-1); % reading line from i position
'R=char(r{1}); %converting cell array to string
u=strfind(R,"); %Finding Location of ' ' in r string
ele=R(1:u(1)-1); %element Name
end
Creating get_mw function
Creating get_mw function
This function takes input as spiecies name and do calculations to evaluate the molecular weight.This function can be understood in following steps:-
a) an array 'a' containing single element name is created and simultaneously another array 'at_w' is created to store element's molecular weight at same index as single element name index.
b) Simpe for loop is run for times equal to length of 'element' array provided at input. Another nested for loop is run for times equal to length of 'a'
c)In nested loop each letter in spiecies name is compared with element name stored in a. 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 is used to find out if spiecies contains any numaric value. Again if statement is used 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.
function [W]=get_mw(element)
a = ['H', 'C' ,'N', '0', 'A'];
at_w = [1 ,12, 14, 16, 40]; W=0;
%Calculating molecular weight of Spieces
for i=1:length(element)
for j=1:length(a)
if (strcmp(element(i),a(j))==1)
W=W+at_w(j);
J=j; %storing the index of the matched element
break; %to terminate For loop of j after strings gets matched
end
end
n=str2double(element(i));
if(n>1)
W=W+(at_w(j).*(n)) -at_w(j);
end
end
end
Creating get_plot function
Creating get_plot function
This function takes inputs as specific heat, enthalpy,entropy , temperature and species name as input.Working of this function can be understood by following statements:-
a)Two strings 'properties of ' and element are merged by strcat command[16] and new directory is created using mkdir having name stored in 'txt' string.
b) similarly txtl ,txt2 and txt2 strings are created and they are used as defining legends title.
c)Then ploting is caried out for specific heat by using simple plot command. To save figure saveas command [17] is used. To decide location of save, directory of program is changed to location stored in 'f. f location string can be set by using fullfile command[18]. After saving, directory location is reset to program directory by using another variable to store intial location of program.
d) Step 'c' is repeated to plot enthalpy and entropy.
function[]=get_plots(cp,H,S,temp,element)
txt=strcat('Properties of ',element); %adding two strings
mkdir(txt); %Making directory of string stored in txt
%Creating Legends title
txtl=strcat(element,' Specific heat vs Temperature');
txt2=strcat(element,' Enthalpy vs Temperature');
txt3=strcat(element,' Entropy vs Temperature');
%Ploting and saving figure(1)
plot(temp,cp);
legend(txtl,'Location','northoutside') % creating legend at particular location
xlabel('Temperature(K)');
ylabel('Specific Heat(Cp) in J/mol*K ');
intial_f=pwd; %storing path of current directory
f=fullfile(txt) %storing path of result directory
cd(f) ; %switching to result directory
saveas(figure(1),[pwd 7specific_heat_vs_Temp_plot.jpg]);
cd(intial_f); %switching back to original directory
%saving figures cd(intial_f); %switching back to original directory
figure(2)
plot(temp,H);
legend(txt2,'Location','northoutside'); % creating legend at particular location
xlabel('Temperature(K)');
ylabel('Enthalpy(h) in 7/mol');
cd(f) ; %switching to result directory
saveas(figure(2),[pwd '/Enthaply_vs_Temp_plot.jpe]); cd(intial_f); %switching back to original directory
figure(3)
plot(temp,S);
legend(txt3,'Location','northoutside'); % creating legend at particular location
xlabel('Temperature(K)');
ylabel('Entropy(s) in 7/mol'K ');
cd(f) ; %switching to result directory saveas(figure(3),[pwd '/Entropy_vs_Temp_plot.jpg']);
cd(intial_f); %switching back to original directory
end
Code
Main code
clc
clear all
close all
format longE; %Set the output format to LongEng.
%opening file thermo30.dat.txt
f1=fopen('thermo30.dat.txt','r');
% Calculating No of lines
for i=1:Inf
line= fgetl(f1);
flag=strcmp(line,'END');
if(flag==1)
no_of_lines=i;
break;
end
end
%finding Global Temperature Range
frewind(f1); % to bring file pointer to intial position of
file fgetl(f1); %reading the first line
data = fgetl(f1); %reading the second line
split = strsplit(data,"); % separating string by space
for i = 1:length(split)
gt(i) = str2double(splitlil); %convrsion of string charaters to numbers
end
global_min = gt(2); %global min temperature
limit = gt(3): %global med temperature
global_high=gt(4); 5global max temperature
%calculating Total number of species present in thermo.data
no_of_species=(no_of_lines-6)/4;
%Data Processing
for k=1:no_of_species
str=k; coff=get_coff(fl,str); %function to get 14 coefficients
h_coff=coff(1:7); % high Temperature coefficients
l_coff=coff(8:14); % low Temperature coefficient
temp=linspace(300,3500,320); % this temperature range is a best fit between local and global temperature
cp=get_cp(h_coff,l_coff,temp,limit); % function to attain cp value for particular Spieces
H=get_h(h_coff,l_coff,temp,limit); % function to attain Enthapy value for particular Spieces
S=get_s(h_coff,l_coff,temp,limit); % function to attain Entropy value for particular Spieces
element=get_element(fl,str); % function to attain spiecies Name value for particular Spieces
molecular_wt(k)=get_mw(element) ; % function to attain spiecies Molecular Weight value for particular Species
get_plots(cp,H,S,temp,element); %function to attain Plots and saving plots to relevant location for particular species
chemical_formula(k)={get_element(fl,str)}; %storing spieces Name in a cell array
end
%generating table to store results of molecular weight
Spieces_Molecular_weight = table(chemical_formula', molecular_wt'); %generating table
Spieces_Molecular_weight.Properties.VariableNames = {'Spieces' Molecular_weight' }; %Naming table variables
writetable(Spieces_Molecular_weight,'Spieces_Molecular_weight.xlsx'); % saving table in xlsx format.
%end
Code function
function[coffl=get_coff(f1,str)
i=4*(str-1)+7; % getting spiecies Number
frewind(f1); % to bring file pointer to intial position of file
d = textscan(f1,'%s',1,'delimiter','n', 'headerlines',i-1); % to read line from i position
D=char(d{1}); %converting cell array to string
E=fgetl(f1); % Reading next line
F=fgetl(f1); % Reading next line
%Finding Location of Exponent in D,E,F string
a=strfind(D,'E');
b=strfind(E,'E');
c=strfind(F,'E');
%getting 14 coefficients
a1=D(1:a(1)+3);
a1=str2double(a1);
a2=D(a(1)+4:a(2)+3);
a2=str2double(a2);
a3=D(a(2)+4:a(3)+3);
a3=str2double(a3);
a4=D(a(3)+4:a(4)+3);
a4=str2double(a4);
a5=D(a(4)+4:a(5)+3);
a5=str2double(a5);
a6=E(1:b(1)+3);
a6=str2double(a6);
a7=E(b(1)+4:b(2)+3);
a7=str2double(a7);
a8=E(b(2)+4:b(3)+3);
a8=str2double(a8);
a9=E(b(3)+4:b(4)+3);
a9=str2double(a9);
a10=E(b(4)+4:b(5)+3);
a10=str2double(a10);
a11=F(1:c(1)+3);
a11=str2double(a11);
al2=F(c(1)+4:c(2)+3);
a12=str2double(a12);
a13=F(c(2)+4:c(3)+3);
a13=str2double(a13);
a14=F(c(3)+4:c(4)+3);
a14=str2double(a14);
coff=[al a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14];
function [Cp]=get_cp(high_temp,low_temp,temp,limit)
R=8.314; %unit j/mol*K
% Calculating Cp for different temperatures
for i=1:length(temp)
terms = high_temp(1);
term2 = high_temp(2).*temp(i);
term3 = high_temp(3).*temp(i)^2;
term4 = high_temp(4).*temp(i)^3;
term5 = high_temp(5).*temp(i)^4;
term6 = low_temp(1);
term7= low_temp(2).*temp(i);
term8 = low_temp(3).*temp(i)^2;
term9 = low_temp(4).*temp(i)^3;
term10 = low_temp(5).*temp(i)^4;
if temp(i)>=limit
Cp(i)= R.*(term1 +term2 +term3 +term4 +term5);
else
Cp(i)=R.*( term6 +term7 +term8 +term9 +term10);
end
end
end
function [H]=get_h(high_temp,low_temp,temp,limit)
R=8.314; %unit 7/mol*K
% Calculating enthalpy for different temperatures
for i=1:length(temp)
term1=high_temp(1).*temp(i);
term2=(high_temp(2).*temp(i)^2)/2;
term3=(high_temp(3).*(temp(i)^3))/3;
term4=(high_temp(4).*(temp(i)^4))/4;
term5=(high_temp(5).*(temp(i)^5))/5;
term6=high_temp(6);
term7=low_temp(1).*temp(i);
term8=(low_temp(2).*temp(i)^2)/2;
term9=(low_temp(3).*(temp(i)^3))/3;
term10=(low_temp(4).*(temp(i)^4))/4;
term11=(low_temp(5).*(temp(i)^5))/5;
term12=low_temp(6);
if temp(i)>=1imit
H(i)=R.*(term1+term2+term3+term4+term5+term6);
else
H(i)=R.*(term7+term8+term9+term10+term11+term12);
end
end
end
function[ele]=get_element(fl,str)
i=4*(str-1)+6; % getting spiecies Number
frewind(f1); % Setting file pointer to intial position of file
r = textscan(f1,'%s',1,'delimiter','n', 'headerlines',i-1); % reading line from i position
R=char(r{1}); %converting cell array to string
u=strfind(R,"); %Finding Location of ' in R string
ele=R(1:u(1)-1); %element Name
end
function [W]=get_mw(element)
a = ['H', 'C' ,'N', '0', 'A'];
at_w = [1 ,12, 14, 16, 40]; W=0;
%Calculating molecular weight of Spieces
for i=1:length(element)
for j=1:length(a) if (strcmp(element(i),a(j))==1)
W=W+at_w(j);
l=j; %storing the index of the matched element
break; %to terminate For loop of j after strings gets matched
end
end
n=str2double(element(i));%finding numbers in spiecies name
if (n>1)
W=W+(at_w(7).*(n)) -at_w(7);
end
end
end
function[]=get_plots(cp,H,S,temp,element)
txt=strcat('Properties of ',element); %adding two strings
mkdir(txt); %Making directory of string stored in txt
%Creating Legends title
txtl=strcat(element,' Specific heat vs Temperature');
txt2=strcat(element,' Enthalpy vs Temperature');
txt3=strcat(element,' Entropy vs Temperature');
%Ploting and saving figure(1)
plot(temp,cp); legend(txt1,'Location','northoutside') % creating legend at particular location
xlabel('Temperature(K)');
ylabel('Specific Heat(Cp) in j/mol*K ');
intial_f=pwd; %storing path of current directory
f=fullfile(txt); %storing path of result directory
cd(f) ; %switching to result directory
saveas(figure(1),[pwd 7specific_heat_vs_Temp_plot.jpg'1); %saving figures
cd(intial_f); %switching back to original directory
figure(2)
plot(temp,H);
legend(txt2,'Location','northoutside'); % creating legend at particular location
xlabel('Temperature(K)');
ylabel('Enthalpy(h) in j/mol');
cd(f) ; %switching to result directory
saveas(figure(2),[pwd '/Enthaply_vs_Temp_plot.jpg']);
cd(intial_f); %switching back to original directory
figure(3)
plot(temp,H);
legend(txt3,'Location','northoutside'); % creating legend at particular location
xlabel('Temperature(K)');
ylabel('Enthalpy(h) in j/mol');
cd(f) ; %switching to result directory
saveas(figure(3),[pwd '/Enthaply_vs_Temp_plot.jpg']);
cd(intial_f); %switching back to original directory
Results
Results for O2
Results for N2
Results for Co2
Table generated
