Parsing NASA thermodynamic data using MATLAB.

Aim: To write a code in Matlab to parse the NASA thermodynamic data file and then calculate the thermodynamic properties of various gas species.

Objective:

• Write a function that extracts the 14 coefficients and calculates the enthalpy, entropy, and specific heats for all the species in the 'THERMO.dat' file. The formulae are shown below.

Here R is the universal gas constant, T is the temperature.

• Calculate the molecular weight of each species and display it in the command window.
• Plot the Cp, Enthalpy, and Entropy for the local temperature range (low temperature: high temperature) specific for each species.
• Save the plots as images with appropriate names and dump them in separate folders for each species with a suitable name. All of these should be generated automatically.
• The main program that you will submit will just take the "THERMO.dat" file as an input and generate Cp, H, and S plots.

Theory:

Data Parsing:

Data parsing is a process of extracting data from a file and use that data for further calculation or study or making some important conclusions and how that information can be used in a computer program.

NASA Thermodynamic Data File:

'THERMO.dat' file is a data file containing 14 coefficients for 53 different elements given by NASA. These coefficients are utilized to calculate the specific heat, enthalpy, and entropy for a particular element. Each element has its own range of local temperature.

MATLAB Program:

Main Program:

% Program to Parse NASA's thermodynamic data
clear all
close all
clc

% Universal Gas constant(R)
R = 8.314;

% loading data file
File = fopen('THERMO.dat','r');

% reads 1st line
fgetl(File);

% collect Global temprature
tline = fgetl(File);
Global_Temp = strsplit(tline,' ');
Gobal_low_temp = str2num(Global_Temp{2});
Gobal_mid_temp = str2num(Global_Temp{3});
Gobal_high_temp = str2num(Global_Temp{4});

% skip comments
fgetl(File);
fgetl(File);
fgetl(File);

for i = 1 : 53  % to iterate through all 53 elements

    % reads species
    line1 = fgetl(File);
    species = get_species(line1);
    
    % calculate molecular weight and print it
    molecular_weight = get_molecular_weight(species);
    fprintf('\nMolecular Weight of %s is : %d',species,molecular_weight)
    
    % Collect local temprature range
    local_temp_range = get_local_temp_range(line1);

    % collect Coefficient data
    a(1:5) = get_coeff(fgetl(File));
    a(6:10) = get_coeff(fgetl(File));
    a(11:14) = get_coeff(fgetl(File));

    % Temprature 
    T = linspace(local_temp_range(1),local_temp_range(3),500);

    % Specific Heat
    Cp = Specific_heat(a, local_temp_range, T, R);

    % Enthalpy
    H = Enthalpy(a, local_temp_range, T, R);

    % Entropy
    S = Entropy(a, local_temp_range, T, R);

    % Plotting and saving in folder

    current_dir = pwd;
    mkdir(species);     % creates new folder with species name
    cd(species);        % open that folder

    % plotting temprature vs specific heat
    figure(1)
    plot(T,Cp,'linewidth',2,'color','b')
    xlabel('Temprature[K]');
    ylabel('Specific Heat[KJ/Kg-K]');
    title(sprintf('Temprature vs Specific Heat plot for %s',species))
    grid on

    % Save jpg of figure(1)
    saveas(figure(1), 'Specific Heat.jpg')

    % plotting temprature vs enthalpy
    figure(2)
    plot(T,H,'linewidth',2,'color','b')
    xlabel('Temprature[K]');
    ylabel('Enthalpy[KJ/Kg]');
    title(sprintf('Temprature vs Enthalpy plot for %s',species))
    grid on

    % Save jpg of figure(2)
    saveas(figure(2), 'Enthalpy.jpg')

    % plotting temprature vs entropy
    figure(3)
    plot(T,S,'linewidth',2,'color','b')
    xlabel('Temprature[K]');
    ylabel('Entropy[KJ/Kg-K]');
    title(sprintf('Temprature vs Entropy plot for %s',species))
    grid on

    % Save jpg of figure(3)
    saveas(figure(3), 'Entropy.jpg')

    cd(current_dir) % goes back to parent folder

end

Function for calculating Molecular Weight:

function molecular_weight = get_molecular_weight(species)

molecular_weight = 0;
Atoms = ['O', 'H', 'N', 'C', 'S', 'A'];
Atomic_weights = [16, 1, 14, 12, 32, 40];

for i = 1 : length(species)
    for j = 1 : length(Atoms)
        if strcmpi(species(i), Atoms(j))
            molecular_weight = molecular_weight + Atomic_weights(j);
            mark = j;
        end
    end
    k = str2num(species(i));
        if k>1
            molecular_weight = molecular_weight + Atomic_weights(mark)*(k-1);
        end
        
    
end

Function for calculating Cp:

function Cp = Specific_heat(a, local_temp_range, T, R)

for i = 1 : length(T)
    if T(i) <= local_temp_range(2)
        Cp(i) = R*(a(8) + a(9)*T(i) + a(10)*T(i)^2 + a(11)*T(i)^3 + a(12)*T(i)^4);
    else
        Cp(i) = R*(a(1) + a(2)*T(i) + a(3)*T(i)^2 + a(4)*T(i)^3 + a(5)*T(i)^4);
    end
end

Function for calculating Enthalpy:

function H = Enthalpy(a, local_temp_range, T, R)

for i = 1 : length(T)
    if T(i) <= local_temp_range(2)
        H(i) = R*T(i)*(a(8) + a(9)*T(i)/2 + a(10)*(T(i)^2)/3 + a(11)*(T(i)^3)/4 + a(12)*(T(i)^4)/5 + a(13)/T(i));
    else
        H(i) = R*T(i)*(a(1) + a(2)*T(i)/2 + a(3)*(T(i)^2)/3 + a(4)*(T(i)^3)/4 + a(5)*(T(i)^4)/5 + a(6)/T(i));
    end
end

Function for calculating Entropy:

function S = Entropy(a, local_temp_range, T, R)

for i = 1 : length(T)
    if T(i) <= local_temp_range(2)
        S(i) = R*(a(8)*log(T(i)) + a(9)*T(i) + a(10)*(T(i)^2)/2 + a(11)*(T(i)^3)/3 + a(12)*(T(i)^4)/4 + a(14));
    else
        S(i) = R*(a(1)*log(T(i)) + a(2)*T(i) + a(3)*(T(i)^2)/2 + a(4)*(T(i)^3)/3 + a(5)*(T(i)^4)/4 + a(7));
    end
end

Function for calculating Species:

function species = get_species(line)

temp = strsplit(line, ' ');
species = (temp{1});

end

Function for calculating Coefficients:

function [a] = get_coeff(line)

E = findstr(line,'E');
a(1) = str2num(line(1:E(1)+3));
for i = 1 : length(E)-1
    temp = str2num(line(E(i)+4:E(i+1)+3));
    a(i+1) = temp;
end

end

Function for calculating Low-Temperature Range:

function [temp] = get_local_temp_range(line)

split = strsplit(line, ' ');
l = length(split);
mid_temp = str2num(split{l-1});
high_temp = str2num(split{l-2});
low_temp = str2num(split{l-3});

temp = [low_temp mid_temp high_temp];
end

Steps:

1) Start the program with 'clear all','close all', and 'clc' commands

2) The only value for input is the value of universal gas constant i.e. 8.314 J/mol-k.

3) NASA's data file 'THERMO.dat' is opened using 'fopen()'.

4) Using 'fgetl()' we move to the next line.

5) For loop is used as the process needs to be repeated for all the 53 elements.

6) Data required like species name, global and local temperatures are collected using commands like 'strsplit()' which is used to split the line into an array, and using indexing of the line we can get the required data.

7) 'str2num()' is used to convert string data into numeric data. It is used in extraction a1 to a14 values and also for finding the number of atoms in the elements like N2. The number of atoms of N is 2.

8) The functions are created for molecular weight and properties (enthalpy, entropy, and specific heat where based on the low-temperature range the coefficients to be used for calculation of the properties are determined).

9) The functions return the values to the main program and the result is printed in the command window.

10) Plots are plotted for temperature vs (enthalpy, entropy, and specific heat) respectively for each and every element and they are saved in separate folders using save as command by changing the directory before saving and switching back to the original directory after saving so that, the process is repeated for all species.

Results:

Command window output:

Folders that were created automatically:

Plots:

O2:

N2:

CO2:

Errors:

No complicated errors occurred while running the MATLAB program.