Calculation of thermodynamic properties for Elements using NASA thermodynamic data

Introduction:

File parsing is a process in which the data’s present is decoded and interpreted with some rules, making them understandable either by the machine or man. Here we are provided with ‘THERMO.dat’ file which contains coefficients and temperature ranges using which we can find out specific heat, enthalpy, and entropy for the species.

The Chemkin 'THERMO.dat' file format is based very closely on the NASA (1971) file format.

Aim:

• To extract the 14 coefficients and calculates the enthalpy, entropy, and specific heats for all the species in the data file.
• Calculate the molecular weight of each species.
• To Plot the Cp, Enthalpy, and Entropy for the local temperature range (low temperature: high temperature) specific for each species.

Governing equation:

The governing equation for finding out the specific heat, enthalpy, and entropy is given below:

Commands used:

1. ‘fopen’- this command is used to open the file that is to be interpreted.
2. The 2nd line in the file represents global temperatures, thus the values are stored for temperature calculations. The 2^nd line is stored into variable ‘Global_temp’
3. The temperature here is inside a string, so we use ‘strsplit’ with a delimiter ‘ ‘. The values are in cell array and thus we take up values using {}. Str2double is used to change the values from string to numerical values.
4. Then skipping the 3 comment lines present in the thermo.dat file, with command 'fgetl(f1)'.
5. Now looping for interpreting the coefficients from the file for each species, every line is extracted as a string using ‘fgetl(f1)’ command and the string ‘E’ positions are found out using ‘strfind’ command.
6. To get a varied temperature range, we use ‘linspace’ command between local high temperature(ht) and local low temperature(lt) with a range of 1000 values.
7. using ‘mkdir’ command, a new directory is created to save the plots for each species in JPG format.
8. After saving the plots again the directory is changed to run the molecular mass script using ‘cd(‘directory’)’ command.

Main parsing program:

clear all 
close all
clc

% universal gas constant
R = 8.314 ; 

%decoding the header infos
f1 = fopen('THERMO.dat','r');
Header= fgetl(f1);

%Reading Global Temperature Line

Global_temp = fgetl(f1);
C = strsplit(Global_temp,' ');
global_low_temp =str2double(C{2});                          %Global Low Temperature 
global_mid_temp =str2double(C{3});                          %Global Mid Temperature 
global_high_temp =str2double(C{4});                          %Global High Temperature 

%%skipping 3 lines
fgetl(f1);
fgetl(f1);
fgetl(f1);

%looping for the next 53 lines

for i = 1:53
line= fgetl(f1);                                        %Reading the line
B=strsplit(line,' ');                                   %Seperating Variable
species=B{1};                                           %Naming First Variable as 'Species'
C = findstr(line,'G');                                  %finding out a string 'G'
local_low_temp = line(C(1)+3:C(1)+6); 
local_high_temp = line(C(1)+13:C(1)+16);
local_mid_temp = line(C(1)+23:C(1)+27);
 
lt = str2double(local_low_temp);                         %Local Low Temperature
ht = str2double(local_high_temp);                       %Local Mid Temperature 
mt = str2double(local_mid_temp);                        %Local High Temperature


%line 1-Reading the coefficient data

co_effs_line1= fgetl(f1);
C = strfind(co_effs_line1,'E');
a1 = co_effs_line1(1:C(1)+3); 
a1 = str2double(a1);                                  %converting str to numerical value
a2 = co_effs_line1(C(1)+4:C(2)+3);
a2 = str2double(a2);
a3 = co_effs_line1(C(2)+4:C(3)+3);
a3 = str2double(a3);
a4 = co_effs_line1(C(3)+4:C(4)+3);
a4 = str2double(a4);
a5 = co_effs_line1(C(4)+4:C(5)+3);
a5 = str2double(a5);

co_effs_line2= fgetl(f1);
D = strfind(co_effs_line2,'E');
a6 = co_effs_line2(1:D(1)+3); 
a6 = str2double(a6); 
a7 = co_effs_line2(D(1)+4:D(2)+3);
a7 = str2double(a7); 
a8 = co_effs_line2(D(2)+4:D(3)+3);
a8 = str2double(a8);
a9 = co_effs_line2(D(3)+4:D(4)+3);
a9 = str2double(a9); 
a10 = co_effs_line2(D(4)+4:D(5)+3);
a10 = str2double(a10); 

co_effs_line3= fgetl(f1);
K = strfind(co_effs_line1,'E');
a11 = co_effs_line3(1:K(1)+3); 
a11 = str2double(a11);
a12 = co_effs_line3(K(1)+4:K(2)+3);
a12 = str2double(a12);
a13 = co_effs_line3(K(2)+4:K(3)+3);
a13 = str2double(a13);
a14 = co_effs_line3(K(3)+4:K(4)+3);
a14 = str2double(a14);


Temp = linspace(lt,ht,1000);

specific_heat = CP(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,ht,lt,mt,R,Temp)  %calling specific heat function

Entrophy = S(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,ht,lt,mt,R,Temp)        %calling entrophy function

Enthalphy = H(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,ht,lt,mt,R,Temp)       %calling enthalpy function



figure(1)
plot(Temp,Entrophy,'linewidth',4,'color','b')
   str = strcat({'Temperature vs Entrophy of'},{' '},{species});         %addding title to the plot
   title(str);
xlabel('Temperature')
ylabel('Entrophy')

figure(2)
plot(Temp,Enthalphy,'linewidth',4,'color','g')
   str = strcat({'Temperature vs Enthalpy of'},{' '},{species});         %addding title to the plot
   title(str);
xlabel('Temperature')
ylabel('Enthalpy')

figure(3)
plot(Temp,specific_heat,'linewidth',4,'color','r')
str = strcat({'Temperature vs Specific heat of'},{' '},{species});%addding title to the plot
title(str);
xlabel('Temperature')
ylabel('Specific heat')

cur_dir=pwd;                                            
mkdir(species);                                         %creating a new directory to save the plots for each species
cd(species);                                            %making the species directory as current one


saveas(1,'Entrophy','jpg');                             %saving the plots in jpg format
saveas(2,'Enthalpy','jpg');
saveas(3,'Specific_heat','jpg');

directory = ('F:/NASA_Datas');              % changing the directory to load the molecular mass script
cd(directory)

MolecularMass = molecular_mass(species);    %calling for the molecular mass function
end

fclose(f1)

function for specific heat:

function [specific_heat] = CP (a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,ht,lt,mt,R,Temp)
if Temp >=1000
    
specific_heat=(a1+a2.*(Temp)+a3.*(Temp.^2)+a4.*(Temp.^3)+a5.*(Temp.^4))*R;

else
specific_heat=(a8+a9.*(Temp)+a10.*(Temp.^2)+a11.*(Temp.^3)+a12.*(Temp.^4))*R;

end

function for Entropy:

function [Entrophy] = S(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,ht,lt,mt,R,Temp)
if Temp>=1000

  Entrophy=(a1.*log(Temp)+a2.*Temp+a3.*(Temp.^2)./2+a4.*(Temp.^3)./3+a5.*(Temp.^4)./4+a7)*R;

else

Entrophy=(a8.*log(Temp)+a9.*Temp+a10.*(Temp.^2)./2+a11.*(Temp.^3)./3+a12.*(Temp.^4)./4+a14)*R;

end

Function for Enthalpy:

function [Enthalphy] = H(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,ht,lt,mt,R,Temp) 
if Temp>=1000
    
Enthalphy=(a1.*Temp+a2.*(Temp.^2)./2+a3.*(Temp.^3)./3+a4.*(Temp.^4)./4+a5.*(Temp.^5)./5+a6)*R;

else
    
Enthalphy=(a8.*Temp+a9.*(Temp.^2)./2+a10.*(Temp.^3)./3+a11.*(Temp.^4)./4+a12.*(Temp.^5)./5+a13)*R;

end

Function for Molecular mass:

function [MM] = molecular_mass(species)
        Atoms=['H','C','N','O','AR'];
        Atomic_Weight=[1 12 14 16 40];
               MM=0;
        for i=1:length(species)
            for j=1:length(Atomic_Weight)
                if strcmp(species(i),Atoms(j))
                    MM=MM+Atomic_Weight(j)
                    u = j;
                end
            end
         % Adding the weights of the repeating atoms
                n=str2double(species(i));
                if n>1
                    MM=MM+Atomic_Weight(u)*(n-1)
                end
        end
      end

Plots of O2:

Plots for N2:

Plots for CO2:

Results:

thus the specific heat, entropy, and enthalpy are calculated for each species.

Temp vs properties graphs is plotted and saved as JPG files for each species in separate directories.

Molecular mass are also calculated.

References:

1.https://skill-lync.com/projects/Parsing-NASA-Thermodynamic-Data-59639

2.https://www.mathworks.com/help/matlab/ref/strsplit.html

3.https://www.mathworks.com/help/matlab/ref/fgetl.html

4.https://www.mathworks.com/matlabcentral/answers/317006-create-new-folders-by-mkdir