MATLAB Program to parse NASA thermodynamic data and calculate thermodynamic properties and molecular weight of various species.

Aim:- To write a program to parse the NASA Thermodynamic data file and then calculate the thermodynamic properties and molecular weight of various gas species.

Governing equations:-

`C_p = R(a_1+a_2*T+a_3*T^2+a_4*T^3+a_5*T^4)`

`H = R*T(a_1+a_2*T/2+a_3*(T^2)/3+a_4*(T^3)/4+a_5*(T^4)/5+(a_6)/T)`

`S = R(a_1*lnT+a_2*T+a_3*(T^2)/2+a_4*(T^3)/3+a_5*(T^4)/4+a_7)`

Where,

`C_p=`Specific Heat

`H=`Enthalpy

`S=`Entropy

`R=`Universal gas constant

`T =`Temperature

`a_1` to `a_7=`Temperature Coefficients

Software used:- MATLAB R2019a

Main code:-

close all 
clear all
clc

%opening the file and reading
f1=fopen('THERMO.dat','r');
c=fgetl(f1);                        %reading the 1st line

%spliting the line and reading the global temperature ranges
A=strsplit(fgetl(f1),' ');          
global_low_temp = str2double(A{2});
global_mid_temp = str2double(A{3});
global_high_temp = str2double(A{4});

%reading the 3 comment lines
for i=1:3
   fgetl(f1);
end

%using while loop to read the file until end
while ~feof(f1)
    %splitting line to get species name 
    line_1=strsplit(fgetl(f1),' ');
    %if the line starts with END the code terminates
    if strcmp(line_1,'END')==1
       break  
    end
    %extracting the species name and local temperature ranges
    species_name=line_1{1};
    l=length(line_1);
    local_low_temp=str2double(line_1{l-3});
    local_mid_temp=str2double(line_1{l-1});
    local_high_temp=str2double(line_1{l-2});
      
    %extracting coefficients by splitiing at the decimal point
         line_2 = fgetl(f1);
         l2=strfind(line_2,'.');
         a1=str2double(line_2(l2(1)-2:l2(1)+12));
         a2=str2double(line_2(l2(2)-2:l2(2)+12));
         a3=str2double(line_2(l2(3)-2:l2(3)+12));
         a4=str2double(line_2(l2(4)-2:l2(4)+12));
         a5=str2double(line_2(l2(5)-2:l2(5)+12));
         
         line_3 = fgetl(f1);
         l3=strfind(line_3,'.');
         a6=str2double(line_3(l3(1)-2:l3(1)+12));
         a7=str2double(line_3(l3(2)-2:l3(2)+12));
         a8=str2double(line_3(l3(3)-2:l3(3)+12));
         a9=str2double(line_3(l3(4)-2:l3(4)+12));
         a10=str2double(line_3(l3(5)-2:l3(5)+12));
         
         line_4 = fgetl(f1);
         l4=strfind(line_4,'.');
         a11=str2double(line_4(l4(1)-2:l4(1)+12));
         a12=str2double(line_4(l4(2)-2:l4(2)+12));
         a13=str2double(line_4(l4(3)-2:l4(3)+12));
         a14=str2double(line_4(l4(4)-2:l4(4)+12));
         
         
    %Range of temperature for which the Cp, H and S is calculated
         T=linspace(local_low_temp,local_high_temp,100);  
         
    %calculating the themal properties
         for i = 1:length(T)
             [Cp(i),H(i),S(i)]=thermal_prop(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,local_mid_temp,T(i));
         end
         
    %making a folder to save the plots with species name
         mkdir(species_name)
         cur_dir = pwd;
         cd(species_name)   
         
    %creating the plots
         figure(1)
         plot(T,Cp,'r','linewidth',2);
         graph_title1=sprintf('Specific heat vs Temperature of %s',species_name);
         title(graph_title1)
         xlabel('Temperature');
         ylabel('Specific heat');
         
         figure(2)
         plot(T,H,'b','linewidth',2);
         graph_title2=sprintf('Enthalpy vs Temperature of %s',species_name);
         title(graph_title2)
         xlabel('Temperature');
         ylabel('Enthalpy');
          
         figure(3)
         plot(T,S,'k','linewidth',2);
         graph_title3=sprintf('Entropy vs Temperature of %s',species_name);
         title(graph_title3)
         xlabel('Temperature');
         ylabel('Entropy');
         
    %saving the plots with following names
         saveas(1,'Cp VS T.png')
         saveas(2,'H VS T.png')
         saveas(3,'S VS T.png')
         
         cd(cur_dir)        %going back to the main directory
         
    %calculating the molecular weight
         molecular_wt=mole_wt(species_name);
    %displaying the molecular wt. in command window
         fprintf(' Molecular weight of %s is %d amu.\n',species_name,molecular_wt);
end

• First the 'THERMO.dat' file was opened using fopen(). And the lines of this file were read and extracted by using the fgetl(). By splitting line at ' ' (space) global low, mid and high temperature were read and stored in global_low_temp, global_mid_temp and global_high_temp.
• Because the next three lines in the file do not have any information as they are comments. There lines were skipped.
• Next the while loop with '~feof' is used to read file until end. Then the line is extracted and splitted at ' '.
• Then by using if statement the extracted line is compared with 'END' (i.e. to check end of file). If so then its exits the while loop as we used break and program stops.
• Then the species name is read and local low, mid & high temperatures are read and stored in resp. variables.
• Then there are 3 lines of coefficients (i.e. 5,5,4) total 14 coefficients. These are also read by fgetl() and later splitted and stored in variables a1,a2,a3,..............a12,a13,a14 by using str2double().
• Then temperature array is defines by linspace() with 100 values between local low and high temperatures. And then Cp, H, S are calculated and stored.
• And the folder is created by the species name and the 3 plots ()are saved in this folder by using saveas() and returned to main directory.
• Then the molecular weight is calculated by giving input as species name.
• The molecular weight along with the corresponding species name is displayed in the command window.

Function codes:-

• Code for calculating thermal properties

function [Cp,H,S]=thermal_prop(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,local_mid_temp,T)

R = 8.314;  %universal gas constant

if (T>=local_mid_temp)
    Cp = R*(a1+a2*T+a3*T^2+a4*T^3+a5*T^4);
    H = R*T*(a1+a2*(T/2)+a3*((T^2)/3)+a4*((T^3)/4)+a5*((T^4)/5)+a6/T);
    S = R*(a1*log(T)+a2*T+a3*((T^2)/2)+a4*((T^3)/3)+a5*((T^4)/4)+a7); 

else
    Cp = R*(a8+a9*T+a10*T^2+a11*T^3+a12*T^4);
    H = R*T*(a8+a9*(T/2)+a10*((T^2)/3)+a11*((T^3)/4)+a12*((T^4)/5)+(a13/T));
    S = R*(a8*log(T)+a9*T+a10*((T^2)/2)+a11*((T^3)/3)+a12*((T^4)/4)+a14); 
    
end

• In this function we find the thermodynamic properties of species.
• First we define universal gas constant 'R'.
• Then we calculate Cp, H, S by using if else statement, if temperature is greater than or equal to mid temperature we use first 7 coefficients else we use second 7 coefficients to calculate.

• Code for calculating molecular weight

function [molecular_weight] = mole_wt(species_name)
 
Species = upper(species_name);
constituents = {'H','C','N','O','AR'};  %atoms 
atomic_weight = [1 12 14 16 40];        %atomic wt. of species
molecular_weight = 0;                   %initially the value is set to 0

%calculating the wt only based on constiuents/atoms
for i = 1:length(Species)
    for j = 1:length(constituents)
        if strcmp(Species(i),constituents(j))
            molecular_weight = molecular_weight + atomic_weight(j);
            u = j;
        end
    end
%calcuting wt. if there are molecules and adding with previous ans    
    n = str2double(Species(i));
    if n > 1
        molecular_weight = molecular_weight + atomic_weight(u)*(n-1);
    end
end

• In this function we calculate the molecular weight of the species.
• We defined the species, constituents and the respective atomic weights. And set the molecular weight to 0.
• Then by using for loop for the length of species name we use another for loop to compare the species name and the constituents and calculate the molecular weight and we also see if there are more atoms of same species (like O2,N2) and calculate the weight according to the no. of atoms.

Output:-

Command window

Window where all the folders are saved

We can see that there are 53 folders (i.e. total no. of species in the given data file) with the 3 plots in it.

Plots:-

`O_2`

`N_2`

`CO_2`

Main code if the user had to input the required species name:-

It is almost same as the main code above but there are slight changes as shown below

• We insert a new line for the input of species name by user using input().
• And before extracting the local temperatures we first compare the input with the species name by using if statement and strcmpi() which is case insensitive.
• If we input any species name which is not in the data file it doesn't give any output.
• And we add another end at the end of code to terminate if statement.
• As shown in comand window i entered "co2" as input and i got the molecular weight and plots as output shown below.
• The part of code which is highlighted is new added for the user input requirement.