PARSING NASA THERMODYNAMIC DATA

PARSING NASA THERMODYNAMIC DATA

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.The code below parses the therm.dat file to make use of the information.

Main function :

• The program uses following commands 'fopen' command opens the thermodynamic data file [therm.dat] and stores it to a fileId called 'f1' . 'fgetl' of f1 command takes the value of first line of the therm.dat file.
• The 2nd line in the file represents global temperatures,thus the values are stored for temperature calculations .
• Since the temperatures are inside the string we use'strsplit' command that splits the string in the specified delimiter ,here ' ' (space) is the delimiter ,after spliting the values are assinged to variables ( global low,mid and high temperatures) ,using 'str2double' command string value is converted to the numeric value.
• A temperature range with 500 values between global high and low temperature is made by using 'linspace' command,the number of values is specified by 'linspace' commad.
• Preallocations are made to store the properties or data of all the elements present in the file.
• A 'for' loop then skips the 3 comment lines present in the therm.dat file,which has command 'fgetl of f1'.
• Command 'cd(pwd)' opens present working directory. 'mkdir' command makes new folder with name specified inside brackets
• A folder named 'Elements' is created in the working directory which is used to store all the related output files of the program .Command 'cd()' makes the new folder as the directory .
• At last a 'for' loop must run 'j' times ,where= j is number of elements.
• After the loop gets executed 'cd ..' is used to get back to the original working directory.

% PARSING NASA thermodynamic data
clear all
close all
clc

%'fopen'command opens the thermodynamic data file [therm.dat] 
f1=fopen('thermo.txt','rt');  

%'fgetl'command takes the value of first line of the therm.dat file
fgetl(f1); 
%the 2nd line in the file represents  global temperatures,thus the values are stored for temperature calculations
TEMP=fgetl(f1);

%'strsplit' command splits the string in the specified delimiter
global_temp=strsplit(TEMP,' ');  %here' '(space) is the delimiter
%after spliting the values are assinged to variables
%'str2double' command is used to convert the string value to the numeric value
global_low_temp=str2double(global_temp{2});
global_mid_temp=str2double(global_temp{3});
global_high_temp=str2double(global_temp{4});

%a temperature range with 500 values between global high and low temperature is made by using 'linspace' command
T=linspace(global_low_temp,global_high_temp,500);

%preallocations are made to store the properties of all the elements present in the file
a=cell(60,20);

%this 'for'loop skips the 3 comment lines present in the therm.dat file
for i=1:3
comment=fgetl(f1);
end

%creates a folder named 'Elements' in the directory which is used to store all the related output files of this program
   cd(pwd)
   Elements='Elements';
   dirname=Elements;
   mkdir(dirname)   %this command makes new folder with name specified
   cd(dirname)      %this command makes the new folder as the directory  

%this 'for'loop must run 'j' times where= j is number of elements
for j=2:54
%Function that calculates all the coefficients of the element  
[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,element]=coefficient(f1);
%Function that calculates Specific_heat[Cp] of the element
Cp=Specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,global_mid_temp);
%Function that calculates Enthalpy[H] of the element
H=Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,global_mid_temp);
%Function that calculates Entrophy[S] of the element
S=Entrophy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,global_mid_temp);
%Function that plots Cp,H and S vs temperature of an element
plots_Cp_H_S(Cp,H,S,T,element);
%Function that calculates Molecular_weight of the element
Molecular_weight=molecular_w(element);
%Function that stores all the data of each element
a=Properties(a,j,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Cp,H,S,Molecular_weight,element);
end
cd ..  %'..'indicates one level up from the current folder

Sub function for calculating Coefficients :

• The function takes all the data from the therm.dat file after the heading line, global temperatures line and the 3 comments , following data after the comments is the data related to the elements, each element has 4 lines of data in same pattern and there are about 53 elements or compounds in the input file all repeated with same pattern.
• For an element 1st line represents the symbol and data related to the element 'fgetl(f1)' command first takes the string of 1st line and stores in variable 'head' after which command 'strok' takes the first charecter of 'head' which is the symbol of the element in the given file and stores in variable 'element'.
• Next 3 lines are the 14 temperature coefficients of the element in which first 7 coefficients are high temperature coeffients and next 7 are low temperature coefficients.Again 'fgetl'command is used 3times to store next 3 lines.
• Next 'strfind'or'findstr'command is used to determine the positions of the letter 'E' in the line. since every coefficient has letter 'E' in the same position ,it is easy to get the values of the coefficients .in a coefficient 'E'is present at one position ,thus value of coefficient is from 1st to 'E'th position followed by next 3 position. The entire string of coefficient is converted to numeric value by 'str2double'or 'str2num' and is stored to variabes a1 to a14 respectively.

%Function calculates all the coefficients of the element  
function [a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,element]=coefficient(f1)

head=fgetl(f1);       %takes the 1st string line
element=strtok(head); %takes the 1st charecter from 1st line
tline1=fgetl(f1);     %takes 2nd line as tline1
tline2=fgetl(f1);     %takes 3rd line as tline2
tline3=fgetl(f1);     %takes 4th line as tline3
b=strfind(tline1,'E'); %finds the letter 'E' position in string

%all coefficients are converted from string to number and 
%are assigned to variabes from a1-a14 
a1=str2double(tline1(1:b(1)+3));
a2=str2double(tline1(b(1)+4:b(2)+3));
a3=str2double(tline1(b(2)+4:b(3)+3));
a4=str2double(tline1(b(3)+4:b(4)+3));
a5=str2double(tline1(b(4)+4:b(5)+3));
a6=str2double(tline2(1:b(1)+3));
a7=str2double(tline2(b(1)+4:b(2)+3));
a8=str2double(tline2(b(2)+4:b(3)+3));
a9=str2double(tline2(b(3)+4:b(4)+3));
a10=str2double(tline2(b(4)+4:b(5)+3));
a11=str2double(tline3(1:b(1)+3));
a12=str2double(tline3(b(1)+4:b(2)+3));
a13=str2double(tline3(b(2)+4:b(3)+3));
a14=str2double(tline3(b(3)+4:b(4)+3));
end

Sub function for calculating Specific Heat :

• The function takes the values of a1-a14 coefficients, 'T' temperature values [500 values] and the global_mid_temperature .
• Universal gas constant 'R' is defined with its values.
• Conditional 'if 'statement is used to execute the following polynomial to find the entrophy in which if the temperature 'T' is more than the global mid temperature then 1st 7 high temperature coefficients are taken to calculate the specific heat else next seven lower temperature coefficients are used. Thus gives the output specific heat with number of values same as temperature values present.

%Function  calculates Specific_heat[Cp] of the element
function Cp=Specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,global_mid_temp)

R=8.314; %Universal gas constant

if T>=global_mid_temp 
%high temperature coefficients used
Cp=R*(a1+(a2.*T)+(a3.*T.^2)+(a4.*T.^3)+(a5.*T.^4)); 
else
%low temperature coefficients used
Cp=R*(a8+(a9.*T)+(a10.*T.^2)+(a11.*T.^3)+(a12.*T.^4));
end
end

Sub function for calculating Enthalpy :

• The function takes the values of a1-a14 coefficients, 'T' temperature values[500 values] and global_mid_temperature .
• Universal gas constant 'R' is defined with its values.
• Conditional 'if 'statement is used to execute the following polynomial to find the entrophy in which if the temperature 'T' is more than the global mid temperature then 1st 7 high temperature coefficients are taken to calculate the enthalpy else next seven lower temperature coefficients are used. Thus gives the output enthalpy with number of values same as temperature values present.

%Function  calculates Enthalpy[H] of the element
function H=Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,global_mid_temp)

R=8.314;   %Universal gas constant

if T>=global_mid_temp
%high temperature coefficients used
H=R*(a1.*T+(a2.*T.^2)/2+(a3.*T.^3)/3+(a4.*T.^4)/4+(a5.*T.^5)/5+a6); 
else
%low temperature coefficients used
H=R*(a8.*T+(a9.*T.^2)/2+(a10.*T.^3)/3+(a11.*T.^4)/4+(a12.*T.^5)/5+a13);
end 
end

Sub function for calculating Entrophy :

• The function takes the values of a1-a14 coefficients, 'T' temperature values[500 values] and global_mid_temperature .
• Universal gas constant 'R' is defined with its values.
• Conditional 'if 'statement is used to execute the following polynomial to find the entrophy in which if the temperature 'T' is more than the global mid temperature then 1st 7 high temperature coefficients are taken to calculate the entrophy else next seven lower temperature coefficients are used. Thus gives the output entrophy with number of values same as temperature values present.

%Function calculates Entrophy[S] of the element
function S=Entrophy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,global_mid_temp)

R=8.314;      % Universal gas constant 

if T>=global_mid_temp
%high temperature coefficients used
S=R*(a1.*log(T)+(a2.*T)+(a3.*T.^2)/2+(a4.*T.^3)/3+(a5.*T.^4)/4+a7);
else
%low temperature coefficients used
S=R*(a8.*log(T)+(a9.*T)+(a10.*T.^2)/2+(a11.*T.^3)/3+(a12.*T.^4)/4+a14);
end
end

Sub function for making plots :

• Function plots Specific Heat,Enthalpy and Entrophy verses temperature of an element respectivey.
• Function takes the inputs Cp, H, S, T and the element ,thus creates the plots as shown.
• Command 'cd' opens the directory 'Element folder created in the working directory' and 'mkdir' creates a new folder with element name in the directory .
• Three figures are used to plot Specific Heat, Enthalpy and Entrophy verses temperature respectively, the plots takes the title with plot name in first line and the symbol of the element in the second line.
• Now 'cd' command is used to open the directory folder with the name of the element that was created previously. It then saves the three plots with particuar plot names.
• The 'cd ..' command is used to close the directory folder with the name of the element and level up from the current folder to the previous directory ie. folder 'Elements'.

%Function  plots Cp,H and S vs temperature of an element
function plots_Cp_H_S(Cp,H,S,T,element)

cd(pwd) %open directory 'Elements'
dir_name=element; 
mkdir(dir_name) %makes new folder(element folder)

% three figures of the plots
figure(1)
plot(T,Cp,'linewidth',2,'color','r')
title({'SPECIFIC HEAT vs TEMPERATURE for';char(element)})
xlabel('Temperature')
ylabel('Specific heat')

figure(2)
plot(T,H,'linewidth',2,'color','r')
title({'ENTHALPY vs TEMPERATURE for';char(element)})
xlabel('Temperature')
ylabel('Enthalpy')

figure(3)
plot(T,S,'linewidth',2,'color','r')
title({'ENTROPHY vs TEMPERATURE for';char(element)})
xlabel('Temperature')
ylabel('Entropy')

%change directory to the element folder
cd(dir_name)
% saving the figures in the element folder
saveas(1,'SPECIFIC_HEAT_vs_TEMPERATURE.png')
saveas(2,'ENTHAPHY_vs_TEMPERATURE.png')
saveas(3,'ENTROPHY_vs_TEMPERATURE.png')

cd .. %'..'indicates one level up from the element folder 
end

Sub function for calculating Molecular weight :

• Function calculates Molecular weight of the element by taking the string 'element' as input.
• A variable called 'Element' is defined that stores the basic elements (the basic elements are those required to make the complete element or the compound as present in therm.dat file).
• Another variable stores the atomic weight of the basic elements in the same order as in the variable 'Element'.
• Two 'for' loops are used which run 'i' and 'j' times ,'for loop of i' is outside the 'for loop of j' thus for one value of 'i' 'j' runs for j times where i=length or number of letters in the string 'element' from input j=length or number of basic elements in the string 'Element'
• Each time the 'for loop of i' runs it takes the letter from the variable 'element' and 'for loop of j' takes the letter from the variable 'Element'.For i=1 it takes first letter from 'element' compares with all the letters in 'Element'.
• Conditional if statement uses command 'strcmp' which compares the two letters ,if they are same then atomic weight of that basic element is taken and gets added to molecular weight which is initialy zero ,it also stores the 'j' value of the letter to 'k'.
• But if the statement is wrong it checks whether its a number, for this it converts the string to number. Again conditional 'if' checks whether the number is greater than one ,if its not a number condition fails,if its true then atomic weight of previous basic element (position is stored in k) gets multipied with (n-1) (because once atomic weight was already added) and adds to molecular weight then 'break' command executes which exits from the for loop 'j'.similarly following letter in 'i' loop is evaluated.
• Entire molecular weight is calculated.

%Function  calculates Molecular_weight of the element
function Molecular_weight=molecular_w(element)

%basic elements that are required to make the complete element or the compound
Element={'H' 'C' 'O' 'N' 'A' 'S'};
%atomic weight of the basic elements in the same order
Atomic_weight=[1 12 16 14 40 32] ;
%initial condition
Molecular_weight=0;

for i=1:length(element)
for j=1:length(Element)
    if strcmp(element(i),Element(j))
    Molecular_weight=Molecular_weight+Atomic_weight(j);
    K=j; %stores the positon
    else 
    n=str2double(element(i));
        if n>1
        Molecular_weight=Molecular_weight+Atomic_weight(K)*(n-1);
        break %exits the inside loop(j loop)
        end
    end
end
end
end

Sub function which stores data in a variable :

• The Function stores all the data of each element, for each element it takes all the previous calculated data as input and stores it in a cell array which is preallocated(a=cell(60,20)) ,thus faster execution of the program is attained due to the preallocation.

%Function  stores all the data of each element
function a=Properties(a,j,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Cp,H,S,Molecular_weight,element)
 % names of the properties calculated
   a{1,1}='Elements\coeff';
   coeffs='a';
   coeff=cell(1,20);
for i=2:15
coeff{i}=[coeffs '' num2str(i-1,'%i')];
a{1,i}=coeff{i};
end
   a{1,16}='Specific heat';
   a{1,17}='Enthalpy';
   a{1,18}='Entrophy';
   a{1,19}='Molecular Weight'; 
% assigning all the calcuated data to the cell
a{j,1}=element;
a{j,2}=a1;
a{j,3}=a2;
a{j,4}=a3;
a{j,5}=a4;
a{j,6}=a5;
a{j,7}=a6;
a{j,8}=a7;
a{j,9}=a8;
a{j,10}=a9;
a{j,11}=a10;
a{j,12}=a11;
a{j,13}=a12;
a{j,14}=a13;
a{j,15}=a14;
a{j,16}=mean(Cp);
a{j,17}=mean(H);
a{j,18}=mean(S);
a{j,19}=Molecular_weight;
end

Some of the plots for elements O2, N2 and CO2 are as shown below :

1.1 'O2' Element plot ENTHAPHY_vs_TEMPERATURE

1.2 'O2' Element plot ENTROPHY_vs_TEMPERATURE

1.3 'O2' Element plot SPECIFIC_HEAT_vs_TEMPERATURE

2.1 'N2' Element plot ENTHAPHY_vs_TEMPERATURE

2.2 'N2' Element plot ENTROPHY_vs_TEMPERATURE

2.3 'N2' Element plot SPECIFIC_HEAT_vs_TEMPERATURE

3.1 'CO2' Element plot ENTHAPHY_vs_TEMPERATURE

3.2 'CO2' Element plot ENTROPHY_vs_TEMPERATURE

3.3 'CO2' Element plot SPECIFIC_HEAT_vs_TEMPERATURE

Screenshot of the folders created for different element names and the image data of the plots stored in folder are shown below

All the data stored in the cell array 'a' (variable) is shown below ,which contains all the coefficients and calculated properties of all the elements present in the input file.

References:

• Youtube
• Matlab help
• Google website
• projects skill-lync