Project -1 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 program that can read the high and low temperature data from the given data file and calculate the Cp, H and S is given below.

The main program is given below,

clear all
close all
clc

%Inputs
R= 8.314;

%Reading the data file
file = fopen('therm.dat','r');

%Reading the first line of the file 
l1 =fgetl(file);

%Getting the global minimum, maximum and mid temperatures
l2=fgetl(file);
temperatures = strsplit(l2);
global_min = str2num(temperatures{2});
global_mid = str2num(temperatures{3});
global_max = str2num(temperatures{4});

%Defining the temperature range 
temp = linspace (global_min,global_max,500);

%Ignoring the comments in the file 
for i=1:3
fgetl(file);
end

%Creating a loop to go through all the 53 species of data in the file 
for i=1:53
    %Reading the first line which includes the information for species
    l=fgetl(file);
    
    %Reading the name of the species
    element_name=strtok(l);
    
    %Reading the three lines which contains 7 high temp coefficients & 7 low temp 
    %coefficients
    l1=fgetl(file);
    l2=fgetl(file);
    l3=fgetl(file);
    
    %Finding the position of E so that we can extract the coefficients
    a=strfind(l1,'E');
    b=strfind(l2,'E');
    c=strfind(l3,'E');

    %extracting high temp coeff
    a1=str2num(l1(1:a(1)+3));
    a2=str2num(l1(a(1)+3:a(2)+3));
    a3=str2num(l1(a(2)+4:a(3)+3));
    a4=str2num(l1(a(3)+4:a(4)+3));
    a5=str2num(l1(a(4)+4:a(5)+3));
    a6=str2num(l2(1:b(1)+3));
    a7=str2num(l2(b(1)+4:b(2)+3));

    %Extracting low temp coeff
    b1=str2num(l2(b(2)+4:b(3)+3));
    b2=str2num(l2(b(3)+4:b(4)+3));
    b3=str2num(l2(b(4)+4:b(5)+3));
    b4=str2num(l3(1:c(1)+3));
    b5=str2num(l3(c(1)+4:c(2)+3));
    b6=str2num(l3(c(2)+4:c(3)+3));
    b7=str2num(l3(c(3)+4:c(4)+3));

    %Calculating specific heat
    c_p = specific_heat(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,R,temp,global_min,global_max,global_mid);
    
    %Calculalting enthalpy
    H = enthalpy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,R,temp,global_min,global_max,global_mid);
    
    %calculating entropy
    S = entropy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,R,temp,global_min,global_max,global_mid);
    
    %Plotting graphs
    plotting(c_p,H,S,temp,element_name)
    
    %Calculating molecular weights
    mol_wt(i) = molecular_weight(element_name)
    
    %Function to write molecular weight to a file 
    weight = mol_wt(i);
    writing_to_file(element_name, weight)
    
end

fclose (file)

The commands used and the code is explained as follows

• Firstly in the program the value of the universal gas constant R is given as 8.314 as an input.
• The given data is then formatted as a therm.dat file.
• Then using the fopen command the file is used to open the file in the read mode specified by 'r' keyword in the fopen command.
• The file is read line by line by using the fgetl command. The first line represents the title if the file so it is ignored.
• The next line has the global maximum, global mid and global minimum values which is read bby using fgetl command and split into the three values by using the strsplit command and converted into numbers.
• Then the temperature array is defined using linspace command as 500 values between the global minimum and global maximmum values.
• The next three lines are comments in the file so they are ignored by using a for loop to iterate for 3 values.
• Then a for loop is created to loop through the 53 species data given in the file and clculate the required values.
• Inside the loop the fgetl command is used to read the next line of the file which represents the name of the species.
• The name of the element is stored in a variable using strtok command.
• Then the next 3 lines of the file is read which contains the high and low temperature data of the particular species.
• The read lines are in the form of a string so they are converted into numbers using the str2num command and stored in specific variables so they can be used for further calculations.
• The low temperatures values are stored in variables a1,a2,a2,a4,a5,a6 and a7 and the high temperature values are stored invariables b1,b2,b3,b4,b5,b6 and b7.
• The next line is the reference to a function called specific heat which calculates and returns the value of specific heat using the high and low temperature values as inputs.
• The next line is the reference to a function called enthalpy which calculates and returns the value of enthalpy using the high and low temperature values as inputs.
• The next line is the reference to a function called entropy which calculates and returns the value of entropy using the high and low temperature values as inputs.
• The next line is the reference to a function that plots graphs using the given inputs.
• The next line is the reference to a function that calculates the molecular weights of the different elements.
• The next line is the reference to a function that writes the molecular weights of diferent elements to a file.
• The NASA thermdynamic data file is closed using the fclose command.

The code for the functions used in the above program is explained below,

1. Function to calculate the Specific heat

%Function to caculate specific heat
function c_p = specific_heat(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,R,temp,global_min,global_max,global_mid)
    if temp<=global_mid
        c_p = R*(a1 + a2.*temp + a3.*(temp.^2) + a4.*(temp.^3) + a5.*(temp.^4))
    else
        c_p = R*(b1 + b2.*temp+b3.*(temp.^2) + b4.*(temp.^3) + b5.*(temp.^4))
end

This function uses the high temperature coefficeints, low temperature coefficients, universal gas constant, global minimum, global mid and global maximum temperature values as the input parameters.

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 select the appropriate equations with high and low temperature coefficients is used to calculate the specific heat value denoted by c_p and return that value.

2. Function to calculate the Enthalpy

%Function to calculate enthalpy
function H = enthalpy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,R,temp,global_min,global_max,global_mid)
    if temp <=global_mid
        H = R.*temp.*(a1 + a2.*(temp./2) + a3.*((temp.^2)./3) + a4.*((temp.^3)./4) + a5.*((temp.^4)./5) + a6./temp)
    else
        H = R.*temp.*(b1 + b2.*(temp./2) + b3.*((temp.^2)./3) + b4.*((temp.^3)./4) + b5.*((temp.^4)./5) + b6./temp)
end

This function uses the high temperature coefficeints, low temperature coefficients, universal gas constant, global minimum, global mid and global maximum temperature values as the input parameters.

Then using the ifelse conditional statement the temp values is compared with the global mid temperature values and the appropriate equation with the proper maximum or minimum temperature coefficients is selected and the enthalpy is calculated and its value represented by H is returned.

3. Function to clculate the Entropy

%Function to calculate entropy
function S =entropy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,R,temp,global_min,global_max,global_mid)
    if temp<=global_mid
        S = R*(a1.*log(temp) + a2.*temp + a3.*(temp.^2)./2 + a4.*(temp.^3)./3 + a5.*(temp.^4)./4 + a7)
    else
        S = R*(b1.*log(temp) + b2.*temp + b3.*(temp.^2)./2 + b4.*(temp.^3)./3 + b5.*(temp.^4)./4 + b7)

end

This function uses the high temperature coefficeints, low temperature coefficients, universal gas constant, global minimum, global mid and global maximum temperature values as the input parameters.

Then using the ifelse conditional statement the temp values is compared with the global mid temperature values and the appropriate equation with the proper maximum or minimum temperature coefficients is selected and the entropy is calculated and its value represented by S is returned.

4. Function to plot the graphs

%Function to plot graphs for all elements
function plotting(c_p,H,S,temp,element_name)
cd(pwd)
dir_name=element_name;
mkdir(dir_name);

%Plotting temperature vs specific heat
figure(1)
plot(temp,c_p,'lineWidth',3)
title(sprintf('Specific Heat vs Temperature \r\n for %s:',element_name))
xlabel('Temperature')
ylabel('Specific Heat')

%Plotting temperature vs specific heat
figure(2)
plot(temp,H,'linewidth',3)
title(sprintf('Enthalpy vs Temperature \r\n for %s:',element_name))
xlabel('Temperature')
ylabel('Enthalpy')

figure(3)
plot(temp,S,'linewidth',3)
title(sprintf('Entropy vs Temperature \r\n for %s:',element_name))
xlabel('Temperature')
ylabel('Entropy')

cd(dir_name)
saveas(1,'c_p_vs_temp.png')
saveas(2,'H_vs_temp.png')
saveas(3,'S_vs_temp.png')

cd ..

end

This function has no output arguments it only uses c_p, H, S, temp, element_name as input arguments and plots the required graphs.

The code of the function is explained below.

• Using the mkdir command a new folder is created in the current folder pwd with the same name as the element name.
• The first figure(1) represents the plots the temperature vs specific heat values.
• The first figure(2) represents the plots the temperature vs enthalpy values.
• The first figure(3) represents the plots the temperature vs entropy values.
• Using the saveas command the three figures are saved in the new folder which is created by the element name.
• Using the cd.. command the folder is changed back to the original folder.
• The function runs for 53 times through the for loop and it creates seperate folders for each species and saves all the plots in the folders.

5. Function to calculate the molecular weight of all elements

%Function for clculating molecular weight
function [w] = molecular_weight(element_name)

s = {'H' 'C' 'N' 'O' 'AR'}; 
%atomic weight of each elements
atomic_wt = [1 12 14 16 40];                                                      
w = 0;                                                                       

for i = 1:length(element_name)  
 for j = 1:length(s)
    if strcmp(element_name(i), s(j))                           
    w= w + atomic_wt(j); 
    %storing the position of matched element
    position = j;                                                         
  end 
end 
   
n = str2num(element_name(i));
if n > 1
                                                      
w = w + atomic_wt(position)*(n-1);                                        
end 
end
end

This function for molecular weight calculations requires element name as input argument and it is saved in the variable 'element' and the output argument is molecular weight.

• A string array is created which contains the name of the atoms.
• A variable atomic_wt is created which contains the atomic weight of respective atoms in the string array.
• Then used strcmp command to compare the name of the element with the string array.
• If there is repeatition of any element, str2num is used to find that and (n-1) is multiplied to atomic_wt of element to calculate molecular weight.

6. Function to write the molecular weight data of all elements to a file

    %Writing molecular weight of elements to a file
    function writing_to_file(element_name,weight)
    dir_name = element_name
    mol_file = fopen('mol_wt.docx' , 'a+')
    fprintf (mol_file, 'Molecular weight of %s is %f \n', dir_name,weight)
    fclose (mol_file)
    end
    

The above function takes the element name and the molecular weights as the data inputs and writes the molecular weight along with the name of the element in a file.

• Firstly a new file is created with name 'mol_wt.docx' and it is opened in append mode denoted by 'a+' so that we can write the molecular weight of all elements in a single file.
• Using the fprintf command the name of the element along with their molecular weight is written to the file.
• The file is then closed.

The file with all the molecular weights of all the elements is given in the below dive link

https://drive.google.com/file/d/1rYEfZCu-UkZN76xKDq9y3XGZ1PqeP4zW/view?usp=sharing

Open the above file in microsoft wordpad to view the file.

Results

The graph plot of one of the elements O2 is given below

The thermodynamic plots for the remaining species are given in the below link,

https://drive.google.com/drive/folders/1ucnqFp-abHm9m8HPzuJlVvod6PpN1XoB?usp=sharing

Rferences:
1. Google

2. Skill-lync