Parsing NASA thermodynamic data

AIM:-

1. Write a function that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the data file. You will be reading this file NASA thermodynamic data.
2. Calculate the molecular weight of each species and display it in the command window.
3. Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature: high temperature) specific for each species.
4. Save the plots as images with appropriate names and dump them in separate folders for each species with suitable name. All these should be generated automatically.
5. The main program that you will submit will just take the "THERMO.dat" file as an input and generate Cp, H and S plots.

THEORY:-

Thermodynamic Data

Thermodynamic data for individual species are necessary for performing equilibrium chemistry and finite-rate chemical kinetics calculations. These data typically provide the thermodynamic properties (enthalpy, entropy, and specific heat) of a species as a function of temperature.

Programs - Statistical thermodynamics

(partition_rotvib.m) Matlab script to compute partition functions and thermodynamic properties (including coefficients for polynomial fits in Cantera cti and NASA formats) from spectroscopic data for heteronuclear diatomic molecules. Spectroscopic data for three diatomic molecules. Nitric oxide (NO); Hydroxyl (OH); Methylidyne (CH).

Procedure:-

1. First using this code(f1=fopen('THERMO.dat','r')) it will open the file and read the data automatically,
2. (fgetl) using this command the it will read the first line and if we type the code again it will read next line.
3. (strsplit(line2,' ');) using this code it split numbers in the line.
4. The 14 coefficients for each species were calculated by using for loop. Each species contains 3 lines in which 14 coefficients are present.
5. The separate folders will automatically created for each species to save the plots for each species in separate folders while executing the program.
6. The value of Universal gas constant R = 8.314J/mol-K was declared.
7. (fclose(f1)) COde indicate the data is closed.

Based upon the 14 co-efficient and the value of R the values of Specific heat, Enthalpy and Entropy were calculated for each species by using the following formulae’s

Where,

T=temperature

R=Universal gas constant

H=enthalpy

S=entropy

Cp=Specific heat

MAIN PROGRAM

clear all
close all
clc

f1=fopen('THERMO.dat','r')

line1=fgetl(f1)
line2=fgetl(f1)

%heading lines
comment1=fgetl(f1)
comment2=fgetl(f1)
comment3=fgetl(f1)

A=strsplit(line2,' ');
global_low_temp=str2num(A{2})
global_medium_temp=str2num(A{3})
global_high_temp=str2num(A{4})

%species information

for i=1:53
  species=fgetl(f1);
  species_name=species(1:6)
  
  line_1=fgetl(f1);
  
  L1=strfind(line_1,'E');
  
  c1=str2num(line_1(1:(L1(1)+3)));
  
  c2=str2num(line_1((L1(1)+4):(L1(2)+3)));
  
  c3 = str2num(line_1((L1(2)+4):(L1(3)+3)));

  c4 = str2num(line_1((L1(3)+4):(L1(4)+3)));

c5 = str2num(line_1((L1(4)+4):(L1(5)+3)));


line_2 = fgetl(f1);

L2 = strfind(line_2,'E');

c6 = str2num(line_2(1:(L2(1)+3)));

c7 = str2num(line_2((L2(1)+4):(L2(2)+3)));

c8 = str2num(line_2((L2(2)+4):(L2(3)+3)));

c9 = str2num(line_2((L2(3)+4):(L2(4)+3)));

c10 = str2num(line_2((L2(4)+4):(L2(5)+3)));


line_3 = fgetl(f1);

L3 = strfind(line_3,'E');

c11 = str2num(line_3(1:(L3(1)+3)));

c12 = str2num(line_3((L3(1)+4):(L3(2)+3)));

c13 = str2num(line_3((L3(2)+4):(L3(3)+3)));

c14 = str2num(line_3((L3(3)+4):(L3(4)+3)));
  
  %creating folder for each species
mkdir(['C:\Users\USER\Desktop\program\NASA parsing\species\',species_name]);

Temperature = linspace(300,5000,100);  

R =8.314;                              %Universal gas constant

%calculation of properties based on the lower and upper temperature
for j = 1:length(Temperature)
    if Temperature(j)<3000
Cp = (c1+c2.*Temperature+c3.*Temperature.^2+c4.*Temperature.^3+c5.*Temperature.^4).*R;
H = (c1+(c2.*Temperature)./2+(c3.*Temperature.^2)/3+(c4.*Temperature.^3)/4+(c5.*Temperature.^4/5)+c6./Temperature).*Temperature.*R;
S = (c1*log(Temperature)+c2.*Temperature+(c3.*Temperature.^2)/2+(c4.*Temperature.^3)/3+(c5.*Temperature.^4)/4+c7).*R;
    else
Cp = (c8+c9.*Temperature+c10.*Temperature.^2+c11.*Temperature.^3+c12.*Temperature.^4).*R;
H = (c8+(c9.*Temperature)/2+(c10.*Temperature.^2)/3+(c11.*Temperature.^3)/4+(c12.*Temperature.^4/5)+c13./Temperature).*Temperature.*R;
S = (c8.*log(Temperature)+c9.*Temperature+(c10.*Temperature.^2)/2+(c11.*Temperature.^3)/3+(c12.*Temperature.^4)/4+c14).*R;
 end 
end

%plotting

%Temperature vs Entropy
figure(1)
plot(Temperature,S,'linewidth',4,'color','b')
title(species_name)
xlabel('Temperature')
ylabel('Entropy')

%Temperature vs Enthalpy
figure(2)
plot(Temperature,H,'linewidth',4,'color','k')
title(species_name)
xlabel('Temperature')
ylabel('Enthalpy')

%Temperature vs Specific heat
figure(3)
plot(Temperature,Cp,'linewidth',4,'color','y')
title(species_name)
xlabel('Temperature')
ylabel('Specific heat')

directory = (['C:\Users\USER\Desktop\program\NASA parsing\species\',species_name]);   %changing the filepath to save the plots
cd(directory)

%saving the plots

saveas(1,'Entropy.jpg');
saveas(2,'Enthalpy.jpg');
saveas(3,'Specific heat.jpg');

cd('C:\Users\USER\Desktop\program\NASA parsing\')       %switching the directory back to read thermodat file for molecular mass calculation

MW = molecular_weight(species_name)

end
fclose(f1)

Function Code

function MW = molecular_weight(species_name)
  
species = upper(species_name); 

Constituents = ['H', 'C', 'N', 'O', 'A'];

Atomic_weight = [1.00794 12.0107 14.0067 15.9994 39.948]; 

MW = 0; 

for i=1:length(species)
    for j=1:length(Constituents)
        if strcmp(species(i),Constituents(j))
           MW = MW + Atomic_weight(j); 
           u = j;
        end
    end
    n = str2num(species(i));
    if n>1
        MW = MW + Atomic_weight(u)*(n-1);
    end
end
end

RESULTS

• plot:- O2

• Temperature vs Entropy

• Temperature vs Enthalpy

• Temperature vs Specific Heat

• PLOT:- N2

• PLOT:- CO2

Command window Moleculer weight values

Created folders

Conclusion

The graphs were plotted of Temperature vs Entropy, Temperature vs Enthalpy and Temperature vs Specific heat and saved in the seperate created folders for each species.

Finally, the values of molecular weights were calculated for each species. For this pupose a different function was created and then called in the parsing code to display the molecular weight in command window.

REFFERENCES

• https://shepherd.caltech.edu/EDL/PublicResources/sdt/thermo.html

• https://in.mathworks.com/help/matlab/ref/fgetl.html