Parsing NASA thermodynamic data using MATLAB

AIM:

Write a code in MATLAB to parse the NASA thermodynamic data file and then calculate thermodynamic properties of various gas species.

OBJECTIVE:

1. Write a function that extracts the 14 coefficient and calculates the enthalpy, entropy and specific heats for all the species in the data file (THERMO.dat).
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. Plots for O₂, N₂ and C0₂ species to be shown.

THEORY:

Data Parsing:

Data parsing is a process of extracting data from a file and use that data for further calculation or study or making some important conclusions and how that information can be used in a computer program. We can split a large data into small data, and we can able to store them so that we can manipulate the data.

NASA has come up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S. They have also documented the coefficients that are required to evaluate these polynomials for 1000+ species. The polynomials are shown below:

Here R is the universal gas constant, T is the temperature.

And a₁, a₂, a₃, a₄, a₅, a₆, and a₇ are the numerical coefficients supplied in NASA thermodynamic file.

In this project we are given a file 'THERMO.dat’. It is a data file containing 14 coefficients for 53 different elements given by NASA. By using data parsing we will extract these coefficients from the file. These coefficients are utilized to calculate the specific heat, enthalpy, and entropy for a particular element. The First 7 coefficients are HIGH-temperature coefficients and the Second 7 coefficients are LOW-temperature coefficients. Also, we will use the LOCAL TEMPERATURE RANGES which is unique for each individual species for calculating, and not the GLOBAL TEMPERATURE RANGE.

PROGRAM:

Programming language used is MATLAB.

Main program: This program extracts the coefficients, calculates the molecular weight, enthalpy, entropy, specific heat for each particular species and plots the graphs. It also creates a folder for each particular species and stores the plots for that species in its folder.

% Program to Parse NASA's thermodynamic data 
clear all
close all
clc 
 
% Universal Gas constant(R)
R = 8.314; 
 
% Loading NASA Thermodynamic data into File
File = fopen('THERMO.dat','r'); 
 
% Reading the first line (Header of the file) 
fgetl(File); 
 
% Collecting the Global Tempratures
tline = fgetl(File);        % Getting the second line which has global temperatures
global_temp = strsplit(tline,' ');      % Splitting the line wherever we find space (' ') so that we can differentiate high , low and medium temperatures
global_low_temp = str2num(global_temp{2});      % Alloting global low temperature
global_mid_temp = str2num(global_temp{3});      % Alloting global medium temperature
global_high_temp = str2num(global_temp{4});     % Alloting global high temperature 
 
% Skipping the next three lines from the file to avoid the unneccessary data 
fgetl(File); 
fgetl(File); 
fgetl(File); 
 
% Creating a file named "Molecular Weights" to write the molecular weight for all species
File2 = fopen('Molecular Weights.txt','w'); 
fprintf(File2,'%8s %8s','Species','Molecular weight');      % Giving an heading before storing the molecular weights
fprintf(File2,'\n');        % Moving to next line 

for i = 1 : 53  % To iterate through all 53 elements 
 
    % Now we start reading the NASA Thermodynamic data from the file 
    
    line1 = fgetl(File);        % Reading next line which contains specie name
    name = strsplit(line1, ' ');        % Splitting the line wherever we find space (' ') so that we get the name
    species = (name{1});        % The name of the specie is present in the first cell after splitting
    
    
    % Calculating the molecular weight And displaying it in command window
    molecular_weight = calculate_molecular_weight(species);     % Calculating molecular weight using the function
    fprintf('Molecular Weight of %s is : %d \n',species,molecular_weight);      % Displaying in the command window
    fprintf(File2,'%8s %8.1f',species,molecular_weight);        % Writing the names and molecular weight into the file
    fprintf(File2,'\n');
    
    % Collect local temprature range
    q1 = findstr(line1,'.');        % Finding "." in the line so that we can get local temperatures  
    local_low_temp = str2num(line1(q1(1)-4:q1(1)));     % Extracting local low temperature 
    local_high_temp = str2num(line1(q1(2)-4:q1(2)));    % Extracting local high temperature    
    local_mid_temp = str2num(line1(q1(3)-4:q1(3)));     % Extracting local mideum temperature 
    
    % Creating row vector containg local temperatures
    local_temp_range = [local_low_temp, local_mid_temp, local_high_temp];    
 
    % Extracting coefficient for the particular species using the function "extract_coeff"
    a(1:5) = extract_coeff(fgetl(File));        % Extracting first 5 coefficients 
    a(6:10) = extract_coeff(fgetl(File));       % Extracting next 5 coefficients
    a(11:14) = extract_coeff(fgetl(File));      % Extracting next 4 coefficients
 
    % Creating Temperature range 
    T = linspace(local_low_temp,local_high_temp,500); 
 
    % Calculating Specific Heat
    Cp = Specific_heat(a, local_temp_range, T, R); 
 
    % Calculating Enthalpy 
    H = Enthalpy(a, local_temp_range, T, R); 
 
    % Calculating Entropy
    S = Entropy(a, local_temp_range, T, R); 
 
    % Plotting and saving in folder 
    
    current_dir = pwd;  % Alloting current working directory to the variable named "current_dir" 
    mkdir(species);     % Creates new folder with species name 
    cd(species);        % Opens that folder       
    
    % Plotting Temprature vs Specific Heat  
    figure(1)   
    plot(T,Cp,'linewidth',2,'color','r')
    xlabel('Temprature[K]');   
    ylabel('Specific Heat[KJ/Kg-K]');  
    title(sprintf('Temprature vs Specific Heat plot for %s',species)) 
    grid on 
 
    % Save jpg of figure(1) 
    saveas(figure(1), 'Specific Heat.jpg') 
 
    % Plotting Temprature vs Enthalpy    
    figure(2)
    plot(T,H,'linewidth',2,'color','k') 
    xlabel('Temprature[K]');   
    ylabel('Enthalpy[KJ/Kg]'); 
    title(sprintf('Temprature vs Enthalpy plot for %s',species))
    grid on 
 
    % Save jpg of figure(2)
    saveas(figure(2), 'Enthalpy.jpg')
    
 
    % Plotting Temprature vs Entropy   
    figure(3)    
    plot(T,S,'linewidth',2,'color','b') 
    xlabel('Temprature[K]');  
    ylabel('Entropy[KJ/Kg-K]');  
    title(sprintf('Temprature vs Entropy plot for %s',species)) 
    grid on 
 
    
    
    % Save jpg of figure(3)    
    saveas(figure(3), 'Entropy.jpg') 
 
    cd(current_dir)     % Going back to parent folder 
 
end 
fclose(File2);       % Closing file named Molecular Weights
fclose(File);        % Closing the THERMO.dat file

Function for calculating molecular weight for each species:

% Function for calculating Molecular weight of the species found

function molecular_weight = calculate_molecular_weight(species) 
 
molecular_weight = 0; 
% Name of atoms 
atom_names = ['O', 'H', 'N', 'C', 'S', 'A'];
% Atomic mass for each atom in order
atomic_mass = [16, 1, 14, 12, 32, 40]; 
 
for i = 1 : length(species) 
    for j = 1 : length(atom_names)   
        if strcmpi(species(i), atom_names(j))        % Using strcmpi for not being case sensitive
            molecular_weight = molecular_weight + atomic_mass(j);     % If a match is found, then adding weight of that atom.
            pos = j;        % Marking the position to know which atom is found  
        end
    end
    n = str2num(species(i)); 
    if n>1      % Checking if n is numeric    
        % Adding the weight if more than one atom of a particular kind is present.
        molecular_weight = molecular_weight + atomic_mass(pos)*(n-1); 
    end
end

Function for extracting coefficients:

% Function for calculating coefficients

function [a] = extract_coeff(line) 
 
occurance = findstr(line,'E');      % Finding the occurance of E in the line
a(1) = str2num(line(1:occurance(1)+3));     % Extracting first coefficient
for i = 1 : length(occurance)-1 
    temp = str2num(line(occurance(i)+4:occurance(i+1)+3));      % Exctracting one coefficient during one iteration 
    a(i+1) = temp;      % Concatinating all coefficients of one line in one vector
end 
 
end

Function for calculating Specific Heat:

% Function for calculating Specific Heat

function Cp = Specific_heat(a, local_temp_range, T, R) 
 
for i = 1 : length(T)
    if T(i) <= local_temp_range(2)
        % Second 7 coefficients are for LOW Temperature 
        Cp(i) = R*(a(8) + a(9)*T(i) + a(10)*T(i)^2 + a(11)*T(i)^3 + a(12)*T(i)^4);
    else
        % First 7 coefficients are for HIGH Temperature 
        Cp(i) = R*(a(1) + a(2)*T(i) + a(3)*T(i)^2 + a(4)*T(i)^3 + a(5)*T(i)^4);
    end
end

end

Function for calculating Enthalpy:

% Function for calculating Enthalpy

function H = Enthalpy(a, local_temp_range, T, R) 
 
for i = 1 : length(T) 
    if T(i) <= local_temp_range(2)  
        % Second 7 coefficients are for LOW Temperature 
        H(i) = R*T(i)*(a(8) + a(9)*T(i)/2 + a(10)*(T(i)^2)/3 + a(11)*(T(i)^3)/4 + a(12)*(T(i)^4)/5 + a(13)/T(i));
    else
        % First 7 coefficients are for HIGH Temperature 
        H(i) = R*T(i)*(a(1) + a(2)*T(i)/2 + a(3)*(T(i)^2)/3 + a(4)*(T(i)^3)/4 + a(5)*(T(i)^4)/5 + a(6)/T(i));
    end
end

end

Function for calculating Entropy:

% Function for calculating Entropy

function S = Entropy(a, local_temp_range, T, R) 
 
for i = 1 : length(T) 
    if T(i) <= local_temp_range(2) 
        % Second 7 coefficients are for LOW Temperature 
        S(i) = R*(a(8)*log(T(i)) + a(9)*T(i) + a(10)*(T(i)^2)/2 + a(11)*(T(i)^3)/3 + a(12)*(T(i)^4)/4 + a(14)); 
    else
        % First 7 coefficients are for HIGH Temperature 
        S(i) = R*(a(1)*log(T(i)) + a(2)*T(i) + a(3)*(T(i)^2)/2 + a(4)*(T(i)^3)/3 + a(5)*(T(i)^4)/4 + a(7)); 
    end
end

end

EXPLANATION OF PROGRAM:

1. The program starts with initialising R with universal gas constant i.e. 8.314 J/mol-k.
2. Then we load the NASA Thermodynamic data from the “THERMO.dat” file to file handle called “File”. We do so by using fopen() function and opening the file in read mode.
3. Then we start with extracting file data, we first take the first line of the file and then we collect the global temperatures which are present in the second line. We do so by Splitting the line wherever we find space (' ') so that we can differentiate high, low and medium global temperatures.
4. We allot the obtained global temperatures to their respective variables. Then we skip the next three lines from the file to avoid the unnecessary data. To skip we use fgetl()
5. Then we create a file named "Molecular Weights" to write the molecular weight for all species into this file and give it a file handle named “File2”.
6. Then we start a for loop to iterate through data of all the 53 elements present in the file.
7. Inside the loop we start with reading next line which contains specie name, Split the line wherever we find space (' ') so that we get the name of the specie, the name of the specie is present in the first cell after splitting and we store it in a variable called “spicies”.
8. Then we calculate the molecular weight of the specie and display it in command window. Also we write these weights into the file Molecular Weights. Calculation is done by a function called “calculate_molecular_weight”.
9. To get the local temperatures we find delimeter “.” in the line so that we can get local temperatures. We store the location of the dots “.” in the variable q1 and with help of these locations we find the local temperatures.
10. Using the positions in q1 and extracting 4 characters before delimeter & 3 characters after delimeter will give us values of tempreature. As extracted data is in string format, str2num() is used to convert them to nos.
11. After extracting we create row vector containing local temperatures.
12. After that we extract coefficient for the particular species using the function "extract_coeff". We pass the next three lines from the file which contain the 14 coefficients, one by one to this function. We store all the extracted variables in the row vector “a”.
13. To calculate thermodynamic properties, we need temperature range. For that we have considered linspace() to create the range in between local & higher local temperature with 500 intervals. This temperature range is stored in the variable “T”.
14. Then we calculate the Specific Heat, Enthalpy and Entropy over the temperature range using functions “Specific_heat”, “Enthalpy” and “Entropy” respectively.
15. After calculating all the values for thermodynamic properties, we plot the graph of these properties against temperature using plot(). The axes of the plots are labelled by using xlabel() & ylabel(). Species name is given as title to plot so as to categorize it as per species using title().
16. After plotting, next objective is to save them to local disk in separate folder as per species name. For this purpose, we change the path of directory of running script to location where we need to create the folders as per species name using cd(). We used mkdir() to create the folder for each species. By this way, we created 53 folders as per species name. After creating folder we have changed the path of directory once again to the folders of species, so that each time we use saveas() to save the plot, that plots will be saved in its respective species folder.
17. After saving all the plots, we change the path of directory once again to original using cd().
18. We then close all the open files. This is done by using fclose().

EXPLANATION OF THE FUNCTIONS USED:

1 . Function for calculating molecular weight for each species:

1. The name of the function is “calculate_molecular_weight”. It takes input as the specie name (species) and gives its molecular weight as output (molecular_weight).
2. We define row vectors “atom_names” and “atomic_mass” which contains the name of the atoms and the mass of the atoms respectively.
3. By using two for loops we check if the input specie name has any of the atom names in it. If a match is found, then we add weight of that atom to the molecular weight. If more than one atom of a particular kind is found its weight is added as many times as its atom is present in the species.

2 . Function for extracting coefficients:

1. The mane of the function is “extract_coeff”. It takes input as one line (line) in the file and gives the output as a row vector (a) containing all the coefficients in that line.
2. It starts with finding the occurrence of “E” in the line, and store the positions of occurrence of “E” in the variable named as “occurance”.
3. With the help of these positions we extract the coefficients and store them in a row vector which is given as output.

3 . Function for calculating Specific Heat:

1. The name of the function is “Specific_heat”. It takes 4 inputs, which are:
   1. A row vector which contains all the coefficients (a).
   2. A row vector which contains the local temperatures (local_temp_range).
   3. Temperature range over which calculation is to be performed (T).
   4. Universal gas constant (R).

And gives an output as a row vector (Cp) which contains the specific heat value for each temperature in the temperature range.

2. It starts with a for loop to iterate over the temperature range for calculation ( for i = 1 : length(T) ).
3. Inside the loop it checks if the particular temperature in the temperature range lies in the LOW temperature (T(i) <= local_temp_range(2)) or HIGH temperature and uses the coefficients accordingly to calculate the Specific Heat using the formula.

4 . Function for calculating Enthalpy:

1. The name of the function is “Enthalpy”. It takes 4 inputs, which are:
   1. A row vector which contains all the coefficients (a).
   2. A row vector which contains the local temperatures (local_temp_range).
   3. Temperature range over which calculation is to be performed (T).
   4. Universal gas constant (R).

And gives an output as a row vector (H) which contains the Enthalpy value for each temperature in the temperature range.

2. It starts with a for loop to iterate over the temperature range for calculation ( for i = 1 : length(T) ).
3. Inside the loop it checks if the particular temperature in the temperature range lies in the LOW temperature (T(i) <= local_temp_range(2)) or HIGH temperature and uses the coefficients accordingly to calculate the Enthalpy using the formula.

5 . Function for calculating Entropy:

1. The name of the function is “Entropy”. It takes 4 inputs, which are:
   1. A row vector which contains all the coefficients (a).
   2. A row vector which contains the local temperatures (local_temp_range).
   3. Temperature range over which calculation is to be performed (T).
   4. Universal gas constant (R).

And gives an output as a row vector (S) which contains the Entropy value for each temperature in the temperature range.

2. It starts with a for loop to iterate over the temperature range for calculation ( for i = 1 : length(T) ).
3. Inside the loop it checks if the particular temperature in the temperature range lies in the LOW temperature (T(i) <= local_temp_range(2)) or HIGH temperature and uses the coefficients accordingly to calculate the Entropy using the formula.

RESULT:

1. Plot of CO₂.

2. Plot of O₂.

3. Plot of N2.

4. Command window showing molecular weights.

5. Screenshot of window where all folders are saved.

6. Molecular weight data, screenshot of Molecular Weights.txt file created during the program.

CONCLUSION:

• Successfully parsed the NASA Thermodynamic data file, using file parsing techniques in MATLAB.
• Calculated the molecular weight for each particular species present in the file and calculated the Specific Heat, Enthalpy and Entropy over the temperature range using the polynomials.
• Plotted graphs for Cp, Enthalpy and Entropy vs the local temperature range (low temperature : high temperature) specific for each species. Saved these plots as images with appropriate names and saved them in separate folders for each species with that species name. All these are generated automatically after running the program.

REFERENCE:

https://in.mathworks.com/matlabcentral/answers/9480-file-handling-between-different-functions

https://in.mathworks.com/help/matlab/ref/str2num.html

https://in.mathworks.com/help/matlab/ref/strsplit.html

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