Parsing NASA Thermodynamic Data
Aim: To write a code in Matlab/Octave to parse the NASA thermodynamic data file and then calculate thermodynamic properties of various gas species.
Objectives:
To write a function that extracts the 14 co-efficients for a given species from the given file.
To write functions that calculate the enthalpy, entropy and specific heats using the NASA polynomials.
To write a function that calculates the molecular weight of each species.
To write a function that plots the Cp, Enthalpy and Entropy of an input species.
Introduction: Parsing is the process of breaking a data block into smaller chunks by following a set of rules, so that it can be more easily interpreted, managed, or transmitted by a computer. NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S using polynomials. They have also documented the co-efficients that are required to evaluate these polynomials for 1000+ species.
Procedure:
The value of the universal gas constant R is given as 8.31 as an input.Then the given THERMO.dat file has to be accessed to obtain the necessary data 'fopen' command is used to open the desired file and the keyword 'r' is used to read the file.This THERMO.dat file is assigned to a variable 'f1'The file is read line by line by using the 'fgetl' command. The first line represents the title of the file so it is ignored.The second line has the global maximum, global mid and global minimum values which is split into the three values by using the strsplit command and converted into numbersAn array of temperature is generated using linspace commandThe next 3 lines are ignored as they contain only commands. A for loop is employed.Again a for loop is created to read the data of all the species. fgetl command inside the for loop makes our task easier Then the next 3 lines of the file is read which contains the high and low temperature data of the particular species.The read lines are in the form of a string so they are converted into numbers using the str2num command and stored in specific variables so they can be used for further calculations.The low temperatures values are stored in variables lc1,lc2,lc3,lc4,lc5,lc6 and lc7 and the high temperature values are stored invariables hc1,hc2,hc3,hc4,hc5,hc6 and hc7.The next part of the program contains a for loop with formulas to calculate Specific Heat, Enthalpy and Entropy for lower temperatures given T<=t.An else command is used to calculate the thermodynamic properties at higher temperature if the temperature fails to satisfy the criteria given abovePlot of Specific Heat vs Temperature is plotted Plot of Enthalpy vs Temperature is plottedPlot of Entropy vs Temperature is plottedUsing the saveas command the three figures are saved in the new folder which is created by the element name.Using the "cd" command the folder is back to the original folder.The function runs for 52 times through the for loop and it creates seperate folders for each species and saves all the plots in the folders.Molecular Weights of different species is calculated.
clear all
close all
clc
%universal gas constant
R=8.31;
f1=fopen('THERMO.dat','r')
line1=fgetl(f1)
%initial lines
for i=1:4
fgetl(f1)
end
T=linspace(300,5000,500);
for i=1:52
line=fgetl(f1)
A=strsplit(line,' ')
species=A(1)
line2=fgetl(f1)
sf=strfind(line2,'E')
cf1=line2(1:sf(1)+3)
cf2=line2(sf(1)+4:sf(2)+3)
cf3=line2(sf(2)+4:sf(3)+3)
cf4=line2(sf(3)+4:sf(4)+3)
cf5=line2(sf(4)+4:sf(5)+3)
line3=fgetl(f1)
cf6=line3(1:sf(1)+3)
cf7=line3(sf(1)+4:sf(2)+3)
cfe1=line3(sf(2)+4:sf(3)+3)
cfe2=line3(sf(3)+4:sf(4)+3)
cfe3=line3(sf(4)+4:sf(5)+3)
line4=fgetl(f1)
cfe4=line4(1:sf(1)+3)
cfe5=line4(sf(1)+4:sf(2)+3)
cfe6=line4(sf(2)+4:sf(3)+3)
cfe7=line4(sf(3)+4:sf(4)+3)
%high temperature coefficients
hc1=str2num(cf1);
hc2=str2num(cf2);
hc3=str2num(cf3);
hc4=str2num(cf4);
hc5=str2num(cf5);
hc6=str2num(cf6);
hc7=str2num(cf7);
lc1=str2num(cfe1);
lc2=str2num(cfe2);
lc3=str2num(cfe3);
lc4=str2num(cfe4);
lc5=str2num(cfe5);
lc6=str2num(cfe6);
lc7=str2num(cfe7);
%low temperature coefficiemts
lc=[lc1 lc2 lc3 lc4 lc5 lc6 lc7];
hc=[hc1 hc2 hc3 hc4 hc5 hc6 hc7];
for j=1:length(T)
t=1000;
if T(j)<=t
Cp(j) = R*(lc(1)+lc(2)*T(j)+lc(3)*(T(j))^2 +lc(4)*(T(j))^3 +lc(5)*(T(j))^4);
H(j) = R*T(j)*(lc(1)+lc(2)*T(j)/2 +lc(3)*((T(j))^2)/3 +lc(4)*((T(j))^3)/4 + lc(5)*((T(5))^4)/5 + lc(6)/T(j));
S(j) = R*(lc(1)*log(T(j))+lc(2)*T(j)+lc(3)*((T(j))^2)/2 +lc(4)*((T(j))^3)/3 + lc(5)*((T(j))^4)/4 + lc(7));
else
Cp(j) = R*(hc(1)+hc(2)*T(j)+hc(3)*(T(j))^2 +hc(4)*(T(j))^3 +hc(5)*(T(j))^4);
H(j) = R*T(j)*(hc(1)+hc(2)*T(j)/2 +hc(3)*((T(j))^2)/3 +hc(4)*((T(j))^3)/4 + hc(5)*((T(5))^4)/5 + hc(6)/T(j));
S(j) = R*(hc(1)*log(T(j))+hc(2)*T(j)+hc(3)*((T(j))^2)/2 +hc(4)*((T(j))^3)/3 + hc(5)*((T(j))^4)/4 + hc(7));
end
end
spec=char(species)
fd = sprintf('%s',spec);
cur_dir=pwd
mkdir(spec);
cd(spec);
%Plot 1, Specific heat (Cp) vs Temperature (T)
figure(1);
plot(T,Cp);
xlabel('Temperature in K');
ylabel('Specific heat in kJ');
title('Specific heat vs Temperature');
%Plot 2, Enthalpy (H) vs Temperature (T)
figure(2);
plot(T,H);
xlabel('Temperature in K');
ylabel('Enthalpy in kJ/(mol)');
title('Enthalpy vs Temperature');
%Plot 3, Entropy (S) vs Temperature (T)
figure(3);
plot(T,S);
xlabel('Temperature in K');
ylabel('Entropy in kJ/(mol K)');
title('Entropy vs Temperature');
% saving the plots that got generated
saveas(figure(1), 'cp.png');
saveas(figure(2), 'h.png');
saveas(figure(3), 's.png');
cd(cur_dir);
end
%Molecular Weight of Particular Species
molecular_weight(i)=Molecular_Weight(species_element);
%Function_for_Molecular_Weight%
function M=Molecular_Weight(species_element)
atoms={'H','C','N','O','AR'};
AtomicWeight=[1 12 14 16 40];
M=0;
for i=1:length(species_element)
for j=1:length(atoms)
if strcmp(species_element(i),atoms(j))
M=M+AtomicWeight(j);
p=j;
end
n=str2num(species_element(i));
if n>1
M=M+(AtomicWeight(p)*(n-1));
end
end
end