MATLAB Program: Parsing NASA thermodynamic data

Aim: Write a program for parsing thermodynamic data.

Objective:

• Write a function to read the data given in the file and extract the coefficients for each species to calculate enthalpy, entropy and specific heats for all the species in the data file.
• Calculate the molecular weight of each species and display it in the command window of MATLAB.
• Plot the specific heat, enthalpy and entropy for the local temperature range specific for each species.
• Write a code to save the plots as images with appropriate names and dump them in separate folders for each species with suitable name.

Theory:

Parsing actually means to split a file or other input into pieces of data that can be easily be stored or manipulated. In file parsing, we split a string into parts then recognize the parts to convert it into something simpler than a string. File parsing is used to read as well as to write data in a file.

NASA's Thermodynamic data file contains coefficients for polynomials that are used to calculate thermodynamic properties for various chemical species. Hence, we will extract these coefficients and perform calculations in MATLAB.

Entropy: Entropy is an extensive property of a thermodynamic system. It is denoted and defined as

`S = (deltaQ)/T`

Its unit is J/K.

Enthalpy: The total enthalpy of a system cannot be measured directly, hence we measure enthalpy change. It is the total energy possesed by the fluid particles in motion and has unit J. Mathematically,

`H = U +PV`

Specific Heat: Amount of heat required to raise unit temperature of a unit mass of a subsatance. Its unit is J/KgK.

Formulae:

`(Cp)/R = a1 + a2 T + a3 T^2 + a4 T^3 + a5T^4`

`H/(RT) = a1 + (a2T)/2 + (a3 T^2)/3 + (a4 T^3)/4 + (a5 T^4)/5 + (a6)/T`

`S/R = a1 lnT + a2 T + (a3 T^2)/2 + (a4 T^3)/3 + (a5 T^4)/4 +a7`

Program:

Main Program

%Parsing of NASA Thermodynamic Data File
clear all
close all
clc

%Reading THERMO data file
f1 = fopen('THERMO.dat','r');
Ist_line = fgetl(f1);

%Reading Global Temperature
Temp = fgetl(f1);

T = strsplit(Temp,' ');

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

%Skipping the comments
for i = 1:3
   fgetl(f1);
end

%Defining Universal Gas constant
R = 8.3145;         %in J/mol-K 

%Reading all species individually
for j = 1:53
    tline = fgetl(f1);
    A = strsplit(tline,' '); 
    species_name = (A{1});
    
    %Reading Global Temperature for each species
    local_low_temp   = str2num(A{end-3});
    local_high_temp  = str2num(A{end-2});
    local_mid_temp   = str2num(A{end-1});
    
    %Reading 1st line of Temperature values
    tline = fgetl(f1);
    a = strfind(tline,'E');
    
    %Converting the string of values to number
    a1 = tline(1:a(1)+3);
    a1 = str2num(a1);
    a2 = tline(a(1)+4:a(2)+3);
    a2 = str2num(a2);
    a3 = tline(a(2)+4:a(3)+3);
    a3 = str2num(a3);   
    a4 = tline(a(3)+4:a(4)+3);
    a4 = str2num(a4);    
    a5 = tline(a(4)+4:a(5)+3);
    a5 = str2num(a5);
    
    %Reading 2nd line of Temperature values
    tline = fgetl(f1);
    a = strfind(tline,'E');
    
    %Converting the string of values to number
    a6  = tline(1:a(1)+3);
    a6  = str2num(a6);
    a7  = tline(a(1)+4:a(2)+3);
    a7  = str2num(a7);
    a8  = tline(a(2)+4:a(3)+3);
    a8  = str2num(a8);    
    a9  = tline(a(3)+4:a(4)+3);
    a9  = str2num(a9);    
    a10 = tline(a(4)+4:a(5)+3);
    a10 = str2num(a10);
    
    %Reading 3rd line of Temperature values
    tline = fgetl(f1);
    a = strfind(tline,'E');
    
    %Converting the string of values to number
    a11 = tline(1:a(1)+3);
    a11 = str2num(a11);
    a12 = tline(a(1)+4:a(2)+3);
    a12 = str2num(a12);
    a13 = tline(a(2)+4:a(3)+3);
    a13 = str2num(a13);    
    a14 = tline(a(3)+4:a(4)+3);
    a14 = str2num(a14);    
    
    %Defining Temperature range for calculations
    T = linspace(local_low_temp,local_high_temp,1000);
    
    %Calculating Specific Heat, Enthalpy, Entropy and Molecular Weight
    Cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,R,T,local_mid_temp);
    
    H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,R,T,local_mid_temp);
    
    S = Entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,R,T,local_mid_temp);
    
    M = molecular_weight(species_name);
    
    %Creating separate Directory to store graphs and values of each species
    cd('D:\\SkillLync\File Parsing\');
    mkdir(species_name);
    cd(species_name);

    %Writing Molecular Weight to a text file
    M_file = fopen('Molecular Weight.txt','w');
    fprintf(M_file,'Molecular Weight of %s is %d',species_name,M);
    fclose(M_file);
    
    %Plotting Temperature VS Specific Heat graph and Saving
    hold off
    figure(1)
    plot(T,Cp,'linewidth',2,'color','r');
    xlabel('Temperature (K)');
    ylabel('Specific Heat (KJ/KgK)');
    title(sprintf('Specific Heat-Temperature Graph of %s',species_name))
    grid on
    saveas(figure(1),'Specific Heat.jpg');
    
    %Plotting Temperature VS Enthalpy graph and Saving
    figure(2)
    plot(T,H,'linewidth',2,'color','g');
    xlabel('Temperature (K)');
    ylabel('Enthalpy (KJ/mol)');
    title(sprintf('Enthalpy-Temperature Graph of %s',species_name))
    grid on
    saveas(figure(2),'Enthalpy.jpg');
    
    %Plotting Temperature VS Entropy graph and Saving
    figure(3)
    plot(T,S,'linewidth',2,'color','b');
    xlabel('Temperature (K)');
    ylabel('Entropy (KJ/K)');
    title(sprintf('Entropy-Temperature Graph of %s',species_name))
    grid on
    saveas(figure(3),'Entropy.jpg');
    
    %Changing Directory to Original Directory and restoring session
    cd('D:\\SkillLync\File Parsing\');
    close all
end
fclose(f1)

Function to calculate molecular weight

%Defining function for Molecular Weight
function M = molecular_weight(species_name)

%Defining atoms and their respective atomic weights
Atoms = ['O','H','C','N','Ar'];
Atomic_weight = [16,1,12,14,40];
M = 0;

%Comparing Atoms to the read data of species name and then calculating
%molecular weight
for i = 1:length(species_name)
    for j = 1:length(Atoms)
        if strcmp(species_name(i),Atoms(j))
           M = M + Atomic_weight(j);
           place = j;
        end
    end
    n = str2num(species_name(i));
    if (n>1)
       M = M + Atomic_weight(place)*(n-1);
    end
end

%Displaying Molecular weight in command window
fprintf('Molecular Weight of %s : ',species_name);
fprintf('%f',M);
disp(' ')
end

Function to calculate Enthalpy

%Defining function for Enthalpy based on the values of Temperature
function H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,R,T,local_mid_temp)

if T >= local_mid_temp
   H = R*T.*(a1+ ((a2*T)/2)+ ((a3*(T).^2)/3)+ ((a4*(T).^3)/4)+ ((a5*(T).^4)/5)+ (a6./T));
else
   H = R*T.*(a8+ ((a9*T)/2)+ ((a10*(T).^2)/3)+ ((a11*(T).^3)/4)+ ((a12*(T).^4)/5)+ (a13./T));
end
end

Function to calculate Entropy

%Defining function for Entropy based on the values of Temperature
function S = Entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,R,T,local_mid_temp)

if T >= local_mid_temp
   S = R*(a1.*log(T)+ a2*T+ ((a3*(T).^2)/2)+ ((a4*(T).^3)/3)+ ((a5*(T).^4)/4)+ a7);
else
   S = R*(a8.*log(T)+ a9*T+ ((a10*(T).^2)/2)+ ((a11*(T).^3)/3)+ ((a12*(T).^4)/4)+ a14);
end
end

Function to calculate Specific Heat

%Defining function for Specific Heat based on the values of Temperature
function Cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,R,T,local_mid_temp)

if T >= local_mid_temp
   Cp = R*(a1+ (a2*T)+ (a3*(T).^2)+ (a4*(T).^3)+ (a5*(T).^4));
else
   Cp = R*(a8+ (a9*T)+ (a10*(T).^2)+ (a11*(T).^3)+ (a12*(T).^4));
end
end

Working of Code:

1. Main program starts with clear all, close all and clc which means clearing workspace, closing all figures and clearing the command window respectively. This is done to ensure that the program is executed smoothly without any garbage values or prompts.
2. The data required for parsing is stored in the file "THERMO.dat". To open this file we use fopen command.
3. Once the file has been opened, we use the following commands. fgetl to read the first line and global temperature is read, strsplit command splits the string at a specified delimiter which in this case is the spacing between the numbers, str2num converts these strings to nmeric value.
4. To skip comments in the file, for loop was used and we defined the universal gas constant.
5. Fo reading the coefficients, we used strfind to find a common term in the string. Here, 'E' is the ccommon term. And these strings are later converted to numeric values.
6. Thus, all the coefficients are extracted and stored to its respective variables. The array for local temperature range is created using linspace. For calculation of specific heat, entropy and enthalpy.
7. Separate functions are created for calculation of specific heat, enthalpy, entropy. In these functions, if else statemant is used to segragate the high temperature and low temperature coefficients.
8. In the for loop, where calculations are being done, directories are made for each species; in these directories the respective molecular weight is stored in a text file, the plots for specific heat, enthalpy and entropy vs temperature are saved with appropriate title and labels. saveas command is used to save the plots.
9. Atlast, figure windows are closed and directory is changed to the parent directory. This ends the for loop.
10. It is a good programming practice to close all files that are still open. In this program, our THERMO.dat file is still open. It is closed by using fclose command. This completes the program for this project.

Outputs:

Output of molecular weight in command window

Sample output plots:

Plot for CO2

NOTE: As there are 53 species in total hence, it is cumbersome to show each plot. Hence, Google drive link has been added in the references section.

Screenshot of the directory where all the folders were created:

Result & Conclusion:

• The data extraction from given file was done succesfully.
• Plotting of the thermodynamic quantities was done successfully.

References:

• MATLAB for Mechanical Engineers - Course provided by Skill Lync
• GNU Octave
• Supporting Material: Thermodynamics - an engineering approach (SIE) - Yunus Cengel, Michael Boles, Mehmet Kanoglu
• Supporting Material: Energy Engineering - VP Katekar Sir

Material: G-Drive link for the code and obtained plots