Project 1- Parsing NASA Thermodynamic Data Using MATLAB

Aim:

To write a Program in MATLAB to parse the NASA thermodynamic data file.

Objective :

To write,

1. a function to extract the 14 coefficients for a given species from the given file.
2. a function to calculate the enthalpy, entropy and specific heats using the NASA polynomials.
3. a function to calculate the molecular weight of each species.
4. a function to plot the Cp, Enthalpy and Entropy of an input species.

Introduction:

File Parsing is the process of analyzing a string of symbols, either in natural language, computer languages or data structures, conforming to the rules of formal grammar.

Within computational linguistics the term is used to refer to the formal analysis by a computer of a sentence or other string of words into its constituents, resulting in a parse tree showing their syntactic relation to each other, which may also contain semantic and other information.

NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S. They have also documented the co-efficients that are required to evaluate these polynomials for 1000+ species.

The thermodynamic properties in terms of the coefficients are given as :

Software used:

MATLAB R2020(online)

Program Code:

Main code:

%Project 1- Parsing NASA Thermodynamic Data Using MATLAB
close all
clear all
clc

%Opening and reading the 'THERMO.dat' file
F = fopen('THERMO.dat','r')

%Extracting the Title from file
line1 = fgetl(F)       

%Extracting the Global temperature values
line2 = fgetl(F) 

%Extracting the Next three comments lines
line3 = fgetl(F)                     
line4 = fgetl(F)   
line5 = fgetl(F)
A = strsplit(line2,' ');  %' ' = delimeteres

global_low_temp =  str2num(A{2})
global_mid_temp =  str2num(A{3})
global_high_temp = str2num(A{4})

%Extracting the species informations for the 53 species in the file
for i=1:53
    
    species= fgetl(F);
     
    species_name = species(1:6)
    Line_1 = fgetl(F);
    
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(F);

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(F);

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 to be stored

mkdir(['/MATLAB Drive/SPECIES_M/',species_name]);

%Temperature range is assumed for calculation
Temperature = linspace(300,5000,100);

%Universal gas constant taken from internet
R =8.314;             

%properties calculation based on lower temperature and upper temperature
for j = 1:length(Temperature)   %for loop 
    if Temperature(j)<3000      %if command
        
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 of if command
end                              %end of for loop

%plotting the values against temperature

%plot for Temperature vs Entropy
figure(1)
plot(Temperature,S,'linewidth',4,'color','r')
title(species_name)
%labelling the axis
xlabel('Temperature')
ylabel('Entropy')

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

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

%changing the filepath to save the plots
directory = (['/MATLAB Drive/SPECIES_M/',species_name]);   
cd(directory)

%saving the plots to folder using saveas command

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

%switching the directory back to read thermodat file for molecular mass calculation
cd('/MATLAB Drive/') 

MM = molecular_mass(species_name)

end
fclose(F)  %closing the opened file

Code for the molecular_mass Function:

% Function to calculate the molecular masses of each individual species

function MM = molecular_mass(species_name)
  species = upper(species_name); 
  Constituents = ['H', 'C', 'N', 'O', 'Ar']; 
  
  Atomic_weight = [1.00794 12.0107 14.0067 15.9994 39.948]; 

MM = 0; 
for i=1:length(species)
    for j=1:length(Constituents)
        if strcmp(species(i),Constituents(j)) %if command
           MM = MM + Atomic_weight(j); 
           
           u = j;   %storing the matched element position to 'u', position = u  
        end         %End of if command
    end
    n = str2num(species(i));
    if n>1
        MM = MM + Atomic_weight(u)*(n-1);
    end
end
end                 %End of the function 

Code Explanation:

1. Initially first enter the basic commands clear all, close all and clc.
2. Then open the given thermodynamic data file using fopen command and read the data in the file and store the data to the variable F.
3. The value of the universal gas constant R is 8.314 .The given data is then formatted as a thermo.dat file.
4. Then the ‘fopen’ command is used to open the file in the read mode (‘r’).
5. Then by using fgetl(), strsplit() and str2num() command the global low ,mid and high temperatures are found in the second line.
6. Then the 14 coefficients for all the species in the given thermodynamic data file are calculated. In data file each species has three lines in which 14 coefficients are present.
7. The file is read line by line by using the ‘fgetl()’ command. All the heading and species related information related in code. The first line represents the title and so it is ignored.
8. The next line has the global high, global mid and global low temperature values which is read by using ‘fgetl()’ command and split into the three values by using the ‘strsplit()’ command and converted into numbers by str2num() command. The next three lines are comments in the file so they are ignored by using fgetl() command
9. Then a ‘for’ loop is created to loop through the 53 species data given in the file and to calculate the required values.
10. Inside the loop, the ‘fgetl’ command is used to read the next line of the file which represents the name of the species. By using the ‘for loop’, the 14 coefficients are calculated for all the species.
11. By using strfind() command, the number of exponential value in first line is found out, by which the number of coefficient in first line can be found. Then the next 3 lines of the file is read which contains the high and low-temperature data of the particular species.
12. The low temperatures values are stored in variables c1, c2, c3 , c4, c5, c6, c7 and the high temperature values are stored invariables c8, c9, c10, c11, c12, c13, c14. By doing the above steps for line two and three, we can calculate the 14 coefficients for all the species.
13. Using the ‘mkdir’ command a new folder is created within the current folder ‘/MATLAB Drive/SPECIES_M/’ (MATLAB online is used by me) with the same name as the species name.
14. Once the coefficients for a species are found, then calculate the entropy, enthalpy and specific heat. First 7 coefficients are low, mid temperature coefficients and next 7 are high temperature coefficients then calculate Cp, H and S.
15. The next is the reference to specific heat(Cp) which calculates and returns the value of specific heat using the high and low-temperature values as inputs. Assuming universal gas constant R=8.314 J/mol-k and taking temperature range from 300 to 5000 and 100 values using linspace() command.

16. Then we use the ‘if else’ conditional statement to check if the temp values are less than or equal to the global mid temperature values and to select the appropriate equations.
17. The next is the reference to enthalpy (H) .This uses the high-temperature coefficients, low-temperature coefficients, universal gas constant, and global minimum, global mid and global maximum temperature values as the input parameters.
    

18. Then it is the reference to a entropy(S) which calculates and returns the value of entropy using the high and low-temperature values as inputs.

19. After the entropy, enthalpy and specific heatfor all the species are calculated, plot the graph between:

i.Temperature vs Entropy

ii.Temperature vs Enthalpy and

iii.Temperature vs Specific heat.

20. Store all the three plots for all the species in the file. The stored file is shown in the above picture

i.The figure(1) represents the plots the temp vs specific heat values.

ii.The figure(2) represents the plots the temp vs enthalpy values.

iii.The figure(3) represents the plots the temp vs entropy values.

21. Using the ‘saveas’ command the three figures are saved in the new folder which is created by the species name.Using the ‘cd’ command the folder is changed back to the original folder.
22. Finally all the three plots for all the species are stored in the file. The function runs for 53 times through the ‘for’ loop and it creates separate folders for each species and saves all the plots in the respective folders of the species.

Function to calculate the molecular weight of all elements: -

1. The molecular weight function requires 'species name' as input argument and it is saved in the variable 'MM' and the output argument is molecular mass.
2. A variable Atomic_weight is created which contains the atomic weight of respective atoms in the string array.
3. Constituents = ['H', 'C', 'N', 'O', 'Ar'];
4. A string array is created which contains the name of the molecules or atoms.
5. Then the ‘strcmp’ command is used to compare the name of the element with the string array.
6. If there is a repetition of an element,then the ‘str2num’ command is used to find it and (n-1) is multiplied to Atomic_weight of the element to calculate molecular weight of the molecule
7. finally the for loop is ended using end command.
8. The program is then executed to get the results.

Output:

For O2:

For N2:

For CO2:

Screenshot of the window where all the folders are saved in the file:

Molecular masses displayed in Command window:

Errors faced:

1)

Forgot to enclose directory with '[ ]' square brackets.

2)

Mentioned the variable name(species_name) incorrectly as "species name" without the underscore.

References:

1. https://in.mathworks.com/help/matlab/ref/fscanf.html
2. https://in.mathworks.com/matlabcentral/answers/196582-best-way-to-parse-text-file

Conclusions:

Finally, the MATLAB program to fulfill the objectives of this project is successfully created and executed to get the required output.