Project 1 - Parsing NASA thermodynamic data

AIM:

To write a MATLAB program that extracts relevant details (14 coefficients and other data) from thermodynamic data submitted by NASA to determine specific heat, enthalpy, entropy and molecular weight of different species.

INTRODUCTION:

NASA came up with polynomials that can be used to evaluate the thermodynamic properties such as specific heat, enthalpy and entropy. They have documented the coefficients that are required to evaluate these polynomials for 1000+ species. The polynomials are as follows:

a1, a2, a3, a4, a5, a6, a7 - High Temperature Coefficients.

a7, a8, a9, a10, a11, a12, a13, a14 - Low Temperature Coefficients.

BASIC DEFINITIONS:

Specific Heat : Specific heat is the amount of heat per unit mass required to raise the temperature of a body by one degree Celsius.

Enthalpy : Enthalpy of a thermodynamic system is the sum of its internal energy and the work required to achieve its pressure and volume.

Entropy : Entropy is the measure of a system's thermal energy per unit temperature that is unavailable for doing some useful work.

Molecular Weight : Molecular weight is the sum of the atomic weight of the atoms constituting a molecule.

PROCEDURE:

1. The file containing NASA thermodynamic data is read using 'fopen' command and assigned to a pointer 'f'.

2. Certain lines of data from the NASA file is skipped after reading since those data are not much useful.

3. The global temperature limits are read and stored in different variables.

4. There are 53 species for which 14 coefficients, species name and local temperature limits are to be extracted from the NASA file. A for loop is used for this purpose.

5. Four separate functions are written for calculating specific heat, enthalpy, entropy and molecular weight. These functions are called into the main program when it is required.

6. The data extracted and calculated are used to plot graphs like specific heat vs temperature, enthalpy vs temperature and entropy vs temperature for all the 53 species. The above mentioned plots for `O_2`, `N_2` and `CO_2` are attached in the report.

7. Finally, the molecular weight of all the 53 species are displayed in the command window.

CODE:

Main Program:

% NASA Thermodynamic Data - File Parsing

clear all
close all
clc

% To Open the External File
f = fopen('THERMO.dat','r');
line_1 = fgetl(f);          %This Line is Skipped
line_2 = fgetl(f);

% To Obtain the Global Temperature Values
global_values = strsplit(line_2,' ');
global_low_temp = str2num(global_values{2});
global_mid_temp = str2num(global_values{3});
global_high_temp = str2num(global_values{4});

line_3 = fgetl(f);          % This line is Skipped
line_4 = fgetl(f);          % This line is Skipped
line_5 = fgetl(f);

% Universal Gas Constant Definition in KJ/K-mol-K
R = 8.314;

for i = 1:53
    T = fgetl(f);
    
    % To extract Species Name
    name_species = T(1:9);
    name_species = deblank(name_species);           % This line is to remove blank spaces
    
    % To extract Local Temperatures
    temp_species = T(49:75);
    local_values = strsplit(temp_species,' ');
    local_low_temp = str2num(local_values{1});
    local_mid_temp = str2num(local_values{3});
    local_high_temp = str2num(local_values{2});
    
    % To Store the Local Temperatures in a Variable
    temp = [local_low_temp local_mid_temp local_high_temp];
    
    temp_range = linspace(local_low_temp,local_high_temp,500);
    
    coeff = fgetl(f);
    coeff1 = findstr(coeff,'E');
    
    % To Extract Coefficients
    a1 = coeff(1:coeff1(1)+3);
    a1 = str2num(a1);
    
    a2 = coeff(coeff1(1)+4:coeff1(2)+3);
    a2 = str2num(a2);
    
    a3 = coeff(coeff1(2)+4:coeff1(3)+3);
    a3 = str2num(a3);
    
    a4 = coeff(coeff1(3)+4:coeff1(4)+3);
    a4 = str2num(a4);
    
    a5 = coeff(coeff1(4)+4:coeff1(5)+3);
    a5 = str2num(a5);
    
    new_line = fgetl(f);
    coeff2 = findstr(new_line,'E');
    
    a6 = new_line(1:coeff2(1)+3);
    a6 = str2num(a6);
    
    a7 = new_line(coeff2(1)+4:coeff2(2)+3);
    a7 = str2num(a7);
    
    a8 = new_line(coeff2(2)+4:coeff2(3)+3);
    a8 = str2num(a8);
    
    a9 = new_line(coeff2(3)+4:coeff2(4)+3);
    a9 = str2num(a9);
    
    a10 = new_line(coeff2(4)+4:coeff2(5)+3);
    a10 = str2num(a10);
    
    new_line1 = fgetl(f);
    coeff3 = findstr(new_line1,'E');
    
    a11 = new_line1(1:coeff3(1)+3);
    a11 = str2num(a11);
    
    a12 = new_line1(coeff3(1)+4:coeff3(2)+3);
    a12 = str2num(a12);
    
    a13 = new_line1(coeff3(2)+4:coeff3(3)+3);
    a13 = str2num(a13);
    
    a14 = new_line1(coeff3(3)+4:coeff3(4)+3);
    a14 = str2num(a14);
    
    
    % To Calculate Specific Heat
    Cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,temp_range,R);
    
    % To Calculate Enthalpy
    H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,temp,temp_range,R);
    
    % To Calculate Entropy
    S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,temp_range,R);
    
    % To Calculate Molecular Weight
    m_weight = weight(name_species);
    
    % To Save Molecular Weight
    species_w = fopen('mass_species.txt','w');
    fprintf(species_w,'Molecular mass of %s is %f',name_species,m_weight);
    fclose(species_w);
    
    % To Plot and Save
    cur_dir = pwd;
    
    % To Make a New Folder for Each Species
    mkdir(name_species);
    cd(name_species);
    
    
    figure(1)
    plot(temp_range,Cp,'LineWidth',2,'Color',[0.2345 0.4567 0.5678]);
    xlabel('Temperature [K]');
    ylabel('Specific Heat [KJ/kg-K]');
    title(name_species);
    
    
    figure(2)
    plot(temp_range,H,'LineWidth',2,'Color',[0.8764 0.9285 0.2378]);
    xlabel('Temperature [K]');
    ylabel('Enthalpy [KJ/kg]');
    title(name_species);
    
    figure(3)
    plot(temp_range,S,'LineWidth',2,'Color',[0.2134 0.9867 0.4898]);
    xlabel('Temperature [K]');
    ylabel('Entropy [KJ/K]');
    title(name_species);
    
    % To Save the Plots for Each Species in Separate Folders
    saveas(1,'Specific Heat.jpg');
    saveas(2,'Enthalpy.jpg');
    saveas(3,'Entropy.jpg');
    
    cd(cur_dir);
    
    % To Print the Molecular Weight of the Species
    sprintf('molecular Weight of %s = %f',name_species,m_weight)
    
end

% To Close the File
fclose(f);

Specific Heat Calculation Function:

% Specific Heat Calculation Function

function [Cp] = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,temp_range,R)

% High Temp Coeff - a1, a2, a3, a4, a5
% Low Temp Coeff - a8, a9, a10, a11, a12

for i = 1:length(temp_range)
    if(temp_range(i)<=temp(2))
        Cp(i) = R*(a8 + (a9*temp_range(i)) + (a10*(temp_range(i).^2)) + (a11*(temp_range(i).^3)) + (a12*(temp_range(i).^4)));
    else
        Cp(i) = R*(a1 + (a2*temp_range(i)) + (a3*(temp_range(i).^2)) + (a4*(temp_range(i).^3)) + (a5*(temp_range(i).^4)));
    end
end
end

Enthalpy Calculation Function:

% Enthalpy Calculation Function

function [H] = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,temp,temp_range,R)

% High Temp Coeff - a1, a2, a3, a4, a5, a6
% Low Temp Coeff - a8, a9, a10, a11, a12, a13

for i = 1:length(temp_range)
    if(temp_range(i)<=temp(2))
        H(i) = R*((a8*temp_range(i)) + (a9*(temp_range(i).^2/2)) + (a10*(temp_range(i).^3/3)) + (a11*(temp_range(i).^4/4)) + (a12*(temp_range(i).^5/5)) + a13);
    else
        H(i) = R*((a1*temp_range(i) + (a2*(temp_range(i).^2/2)) + (a3*(temp_range(i).^3/3)) + (a4*(temp_range(i).^4/4)) + (a5*(temp_range(i).^5/5)) + a6));
    end
end
end

Entropy Calculation Function:

% Entropy Calculation Function

function [S] = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,temp_range,R)

% High Temp Coeff - a1, a2, a3, a4, a5, a7
% Low Temp Coeff - a8, a9, a10, a11, a12, a14

for i = 1:length(temp_range)
    if(temp_range(i)<=temp(2))
        S(i) = R.*((a8*log(temp_range(i))) + (a9*(temp_range(i))) + (a10*(temp_range(i).^2/2)) + (a11*(temp_range(i).^3/3)) + (a12*(temp_range(i).^4/4)) + a14);
    else
        S(i) = R.*((a1*log(temp_range(i))) + (a2*(temp_range(i))) + (a3*(temp_range(i).^2/2)) + (a4*(temp_range(i).^3/3)) + (a5*(temp_range(i).^4/4)) + a7);
    end
end
end

Molecular Weight Calculation Function:

% Molecular Weight Calculation Function

function [m_weight] = weight(name_species)

element = {'H','O','N','Ar','C','S'};
atomic_mass = [1.008 15.999 14.007 39.948 12.011 32.065];
m_weight = 0;

for i = 1 : length(name_species)
    for j = 1: length(element)
        if strcmp(name_species(i),element(j))
            m_weight = m_weight + atomic_mass(j);
            a = j;
        end
        
    end
    
    b = str2num(name_species(i));
    if b>1
        m_weight = m_weight + atomic_mass(a)*(b-1);
    end
end

end

RESULTS:

OXYGEN:

NITROGEN:

CARBON DIOXIDE:

Command Window Output:

In the same way, the molecular weight is displayed on the command window for the remaining species.

Folders containing the Plots for All Species:

ERRORS:

I got confused with the number of brackets '()' put in the calculation functions. For every calculation function, there was an error of same kind as shown below. By checking for the missing brackets and extra brackets, this problem was solved. There were little typos here and there since this was a long program. The errors were rectified just by visual inspection.

Other than these errors, no major errors were encountered in the coding process.

CONCLUSION:

Hence the file containig NASA thermodynamic data is used to perform the required task and the results are verified.

REFERENCES:

1. https://in.mathworks.com/help/matlab/ref/fscanf.html

2. https://in.mathworks.com/help/matlab/ref/deblank.html

3. https://in.mathworks.com/help/matlab/ref/cd.html

4. https://in.mathworks.com/help/matlab/ref/mkdir.html?searchHighlight=mkdir&s_tid=doc_srchtitle