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 numbers
• An array of temperature is generated using linspace command
• The 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 above
• Plot of Specific Heat vs Temperature is plotted
• Plot of Enthalpy vs Temperature is plotted
• Plot of Entropy vs Temperature is plotted
• Using 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 and
  

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

The plots obtained for O2 is as follows: