WRITING A CODE TO PARSING THE NASA THERMODYNAMIC DATA FILE BY USING MATLAB

% Writing a program code for to parising the nasa thermodynamic file by using Matlab

To parsing the NASA thermodynamic data file and then calculate the thermodynamic Properties of various gas species.

The term \'File Parising\'means how the information can read from a particular text file or data file.In matlab extracts the 14 co-efficients Here we are going to see NASA Thermodynamic file which has 53 Constant monoautomic gas species data with low & High temperatures that can be used to evaluate thermodynamic properties such as Specific heat at constant pressure (Cp), Enthalpy(H), and Entropy(S).The formulae are shown below

%STEP:1
%input
%Universal gas constant
R = 8.314; 
%STEP:2
%Reading the data file
file = fopen(\'thermo.dat\',\'r\');
%STEP:3
%Reading the first line of the file
l1 = fgetl(file);
%Getting the global minimum, maximum and mid temperatures
l2 = fgetl(file);
%STEP:4
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,global_mid);
%STEP:5
  %Ignoring the comments in the file 
for i = 1:3;
    fgetl(file);
end
%STEP:6
%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 informatiom for species 
    l0=fgetl(file);
    %STEP:6
    %Reading the name of the species
    species = strtok(l0);
    %Extracting loop
    a = strfind(l0, \'G\');
    global_min = l0(a+1:end-1);
    global_max = str2num(global_min);
    %STEP:7
    %Reading the three lines which contains 7 high temp coefficient & 7 low temperature coefficients 
    l1=fgetl(file);
    l2=fgetl(file);
    l3=fgetl(file);
    %STEP:8
    %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)+4: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));    
    %STEP:9
    %calculating specific heat
    c_p = specific_heat(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,temp,R,global_min,global_mid,global_max);
    %calculating enthalpy
    H = enthalpy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,temp,R,global_min,global_mid,global_max);
    %calculating entropy
    S = entropy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,temp,R,global_min,global_mid,global_max);
    %calcualtion of cp,H and S for each Species
    for i = 1:length(temp)
        % if temp < global_mid_temp
c_p = R*(a1 + a2.*temp + a3.*(temp.^2) + a4.*(temp.^3) + a5.*(temp.^4)+ a6.*(temp.^5)+ a7.*(temp.^6));
S = R*(a1.*log(temp)+ a2.*(temp)+ a3.*(temp.^2)./2 + a4.*(temp.^3)./3 + a5.*(temp.^4)./4 + a6.*(temp.^5)./5 + a7);
H = R*(a1.*(temp)+[a2.*(temp.^2)/.2]+ [a3.*(temp.^3)./3] + [a4.*(temp.^4)./4] + [a5.*(temp.^5)./5] + [a6.*(temp.^6)./6] + a7);
 %else
   c_p = R*(b1 + b2.*temp+b3.*(temp.^2) + b4.*(temp.^3) + b5.*(temp.^4)+ b6.*(temp.^5)+ b7);
   S = R*(b1.*log(temp)+ b2.*(temp)+ b3.*(temp.^2)./2 + b4.*(temp.^3)./3 + b5.*(temp.^4)./4 + b6.*(temp.^5)./5 + b7);
   H = R*(b1.*(temp)+[b2.*(temp.^2)/.2]+ [b3.*(temp.^3)./3] + [b4.*(temp.^4)./4] + [b5.*(temp.^5)./5] + [b6.*(temp.^6)./6] + b7);
    end  
   [W] = molecularweight(species);
   mkdir(species)
   cd(species)
   %STEP:10
   %Plotting 
    % Specific Heat Vs Temperature  
  figure (1)
    plot(temp,c_p,\'linewidth\',3,\'color\',\'r\')
    title(sprintf(\'Specific Heat Vs Temperature plot of %s\', species));
    xlabel(\'Temperature[K]\');
    ylabel(\'specific heat Cp [J/K]\');
    legend(\'specific heat\');
    grid on 
    %Plotting 
    %Enthalpy Vs Temperature
    figure (2)
    plot(temp,H,\'linewidth\',3,\'color\',\'b\')
    title(sprintf(\'Enthalpy Vs Temperature plot of %s\', species));
    xlabel(\'Temperature [K]\');
    ylabel(\'Enthalpy H [J/Kg]\');
    legend(\'Enthalpy\');
    grid on
    %Plotting 
    %Entropy Vs Temperature 
    figure (3)
    plot(temp,S,\'linewidth\',3,\'color\',\'g\')
    title(sprintf(\'Entropy Vs Temperature plot of %s\', species));
    xlabel(\'Temperature [K]\');
    ylabel(\'Entropy S [J/K]\');
    legend(\'Entropy\');
    grid on 
 %graph
saveas(1,\'Specific_heat.png\')
saveas(2,\'Enthalpy.png\')
saveas(3,\'Entropy.png\')
   cd .. 
end

STEP:1

In thermodynamic

Here the input gas constant for R = 8.314

R stands for Universal gas constant

STEP:2

%Reading the data file

Now we are going to read the file name of THERMO.dat which is stored in our program folder.We are using command called as fgetl() or fopen() to read and get the data of the file. we have to ensure THERMO.dat file & program file should be in same folder.

STEP:3

Here we are going to read the first two lines with help of command fgetl.

it consist of to read the \'THERMO\' in first line and to read our temperature limits in second line.

STEP:4

Defining the temperature ranges

here using command called strsplit() is used to read and plot the values of temperature in the second line of dat file.

It stored in varaiable temperature and splitted in three groups named as global min,global mid,global max temperatures.

here we are bringing one more command called \'str2num\'.It is used to convert the values of the temperature from a string to a number.

STEP:5

Ignoring the comments in the Thermo file

Here we want to skip the next three lines which is unwanted for us while writing code in matalb.

For Skipping lines we are using \'FOR\' loop to skip the lines

STEP:6
Here we are creating a For loop to go through all the 53 species of data in the file

STEP:6
Here we are Reading the name of the species and \'strtok\' to return the species names in text string.

we have to save the species name in the cell array in which the species which are returned as text string is stored.
STEP:7
here Read the three lines which contains 7 high temp coefficient & 7 low temperature coefficients with using command called fgetl().

STEP:8
here we are Finding the position of E so that we can extract the coefficients with help of \'findstr()\'command.

\'findstr()\'command is used to find the\'Exponential\' which is common in all the species. we are using strfind(l,\'E\');
we are extracting both high & Low temp coeff
so we try to extract all the \'E\' in each linie. it can be done by (a1=l(1:12+3)) it found that it is in the 12th position of the string.Again \'str2num\'command is used the purpose of this command explained in above .Totally we have to do for all 14 constant values of each monoautomic gas.

For this monoautomic gas we divided in to two coefficients one is high temperature coefficients (a1,a2,a3,a4,a5,a6,a7) and other one is low temperature coefficients (b1,b2,b3,b4,b5,b6,b7).
%STEP:9

In this we are going to calculating the specific heat,enthalpy,entropy and molecularweight and saving all species in one folder to that we are using command called\'mkdir\' &\'cd\'.

\'mkdir()\' is used to create new folder to save the plots and \'cd()\' command is used to change the current folder to create new folder.

%STEP:10

In this we are going to Plotting the different results along with species of monoautomic gas with specific heat, Enthalpy and Entropy Vs Temperature and using \'saveas()\'command to save the plots in form of png.

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

1.Function to calculate 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)+ a6.*(temp.^5)+ a7.*(temp.^6))
else
    c_p = R*(b1 + b2.*temp+b3.*(temp.^2) + b4.*(temp.^3) + b5.*(temp.^4)+ b6.*(temp.^5)+ b7.*(temp.^6))
end

using this function denotes the high temperature coefficients, low temperature coefficients, universal gas constant ,global maximum,global minimum and global mid temperature values are the inout parameters for this function.

Then we use the if else conditional statement to check the if temp values are less then to the mid temp values and the aprropriate equations with high and low temperature coefficients is used to calculate the specific value denoted by c_P and return that value.

2.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*(a1.*(temp)+[a2.*(temp.^2)/.2]+ [a3.*(temp.^3)./3] + [a4.*(temp.^4)./4] + [a5.*(temp.^5)./5] + [a6.*(temp.^6)./6] + a7);
else
    H = R*(b1.*(temp)+[b2.*(temp.^2)/.2]+ [b3.*(temp.^3)./3] + [b4.*(temp.^4)./4] + [b5.*(temp.^5)./5] + [b6.*(temp.^6)./6] + b7);
end

using this function denotes the high temperature coefficients, low temperature coefficients, universal gas constant ,global maximum,global minimum and global mid temperature values are the inout parameters for this function.

Then we use the if else conditional statement to check the if temp values are less then to the mid temp values and the aprropriate equations with high and low temperature coefficients is used to calculate the Enthalpy denoted by H and return that value.

3.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 + a6.*(temp.^5)./5 + a7);
else
    S = R*(b1.*log(temp)+ b2.*(temp)+ b3.*(temp.^2)./2 + b4.*(temp.^3)./3 + b5.*(temp.^4)./4 + b6.*(temp.^5)./5 + b7);
end

using this function denotes the high temperature coefficients, low temperature coefficients, universal gas constant ,global maximum,global minimum and global mid temperature values are the inout parameters for this function.

Then we use the if else conditional statement to check the if temp values are less then to the mid temp values and the aprropriate equations with high and low temperature coefficients is used to calculate the Entropy denoted by E and return that value.

4.Function to calculate Molecular weight:

function W = molecularweight(species)
s = [\'H\' \'C\' \'N\' \'O\' \'S\' \'A\'];
%atomic weight of each species
atomic_weight = [1 12 14 16 32 40];
W = 0;
for i = 1:length(species)
    for j = 1:length(s)
        if strcmpi(species(i),s(j))
            W = W + atomic_weight(j);
            position = j;
        end
    end
    n = str2num(species(i));
    if n>1
        W = W + atomic_weight(position)*(n-1)
          
    end
end
end

In this function for molecular weight calculations requires species as an input argument and it is saved in the variable \'species\' and the output argument is molecular weight.

• A string array is created which contains the name of the atoms variable \'s\'.
• weight of the each species variable atomic_weight created the values of atoms.
• using for loop command and strcmp command used to compare the name of the species with the string array.
• str2num is used to find that and (n-1) is multiplied to atomic_weight species to calculate molecular weight.

Results

1.Screenshot of the directory where all the plots were saved automatically

2. Plot for O2 gas properties:

a) Specific heat(Cp)

b) Enthalpy(H)

c)Entropy(S)

3.Plots for C02 gas properties:

a)Specific Heat(Cp)

b)Enthalphy(H)

c)Entropy(S)

4.Plots for N2 gas properties:

a)Specific heat(Cp)

b)Enthalphy(H)

c)Entropy(S)