Challenge

Parsing NASA thermodynamic data using MATLAB

VVinayak Jauhari4 Jun 2020

Aim – Parsing NASA’s thermodynamic data using MATLAB.

Objective –

To parse NASA’s thermodynamic data file.
To write a function that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the data file.
To calculate the molecular weight of each species and to display it in the command window.
To plot the graphs of specific heat, enthalpy, entropy for the local temperature range.

Theory –

NASA came up with polynomials that can be used to evaluate thermodynamic properties such as specific heat (Cp), enthalpy (H) and entropy (S). They have also documented the co-efficients that are required to evaluate these polynomials for 1000+ species.

We need to write a function that extracts the 14 co-efficients and calculates the enthalpy (H), entropy (S) and specific heats (Cp) for all the species in the data file. The formulae are shown below-

Here R is the universal gas constant (8.314 J/mol-K), T is the temperature in Kelvin (K).

In the given formula, there are 7 numerical co-efficients (a1, a2, a3, a4, a5, a6 and a7). But in the NASA’s thermodynamic file there are two local temperature ranges – low temperature range (300 - 1000) and high temperature range (1000 - 5000). The co-efficients from a1 to a7 are numerical co-efficients of ‘high temperature range’ and the co-efficients from a8 to a14 are numerical co-efficients of ‘low temperature range’.

File Parsing – Parsing means to convert data in a certain format into a more usable format. Parsing refers to the way that human beings analyze a sentence or phrase (in spoken language or text) "in terms of grammatical constituents, identifying the parts of speech, syntactic relations, etc.

Why do we even need to parse files?

In an imaginary world where all data existed in the same format, one could expect all programs to input and output that data. There would be no need to parse files. However, we live in a world where there is a wide variety of data formats. Some data formats are better suited to different applications. An individual program can only be expected to cater for a selection of these data formats. So, inevitably there is a need to convert data from one format to another for consumption by different programs. Sometimes data is not even in a standard format which makes things a little harder.

Code –

Code for specific heat (Cp) function –

% Specific Heat function 

function c_p = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,R,T,global_mid_temp)

if T>global_mid_temp
    c_p = R.*(a1 + a2.*T + a3.*(T.^2) + a4.*(T.^3) + a5.*(T.^4));

else
    c_p = R.*(a8 + a9.*T + a10.*(T.^2) + a11.*(T.^3) + a12.*(T.^4));
        
end

Code for Enthalpy (H) function –

% Enthalpy function 

function H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,R,T,global_mid_temp)

if T > global_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

Code for Entropy (S) function –

% Entropy function 

function S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,R,T,global_mid_temp)

if T > global_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)./5 + a14);
  
end

Code for plotting these thermodynamic properties –

% function for plotting the properties

function plotting(Cp,H,S,T,Element)
%Saving as a folder for different elements
cur_dir = pwd;
mkdir(Element);
cd(Element);

% Plot 1, Specific heat (Cp) vs Temperature (T)
figure(1)
plot(T,Cp,'linewidth',3,'color','b');
xlabel('Temperature [K]');
ylabel('Specific Heat [KJ/(kg-K)]');
title(sprintf('Specific Heat vs Temperature Plot of %s',Element));
saveas(figure(1),'Cp_vs_T.jpg');

% Plot 2, Enthalpy (H) vs Temperature (T)
figure(2)
plot(T,H,'linewidth',3,'color','g');
xlabel('Temperature [K]');
ylabel('Enthalpy [kJ/mol]');
title(sprintf('Enthalpy vs Temperature Plot of %s',Element));
saveas(figure(2),'H_vs_T.jpg');

% Plot 3, Entropy (S) vs Temperature Range (T)
figure(3)
plot(T,S,'linewidth',3,'color','r');
xlabel('Temperature [K]');
ylabel('Entropy in [J/(mol-K)]');
title(sprintf('Entropy vs Temperature Plot of %s',Element));
saveas(figure(3),'S_vs_T.jpg');

cd(cur_dir);

end

Code for molecular weight of elements –

% function for molecular weight

function [MW] = molecular_weight(Element)

elements = ['H' 'C' 'N' 'O' 'Ar'];

% atomic weight of elements
atomic_wt = [1 12 14 16 40];

MW = 0;

for i = 1:length(Element)
    for j = 1:length(elements)
        if strcmp(Element(i), elements(j))
            MW = MW + atomic_wt(j);
            
            % for position of matched elements
            p = j;
            
        end
        
    end
    n = str2num(Element(i));
    if n > 1
        MW = MW + atomic_wt(p)*(n-1);
    end
end

end

Code to write a file for molecular weight of elements –

% function to write molecular weight of elements

function writefile(Elements,molecular_weight)
Molecular_list = fopen('Molecular_Weight.txt','a+');
fprintf( Molecular_list,'Molecular Weight of %s is %d n',Elements,molecular_weight);
fclose(Molecular_list);
end

Code to file parse NASA’s thermodynamic data –

clear all
close all
clc

% Universal gas constant in [J/mol.K]
R = 8.314;

% Opening the data file
f1 = fopen('THERMO.dat','r');      % r is for reading

% For first line
L1 = fgetl(f1);

% For second line
% To get the global min, mid, max temperatures
L2 = fgetl(f1);
A = strsplit(L2);
global_min_temp = str2num(A{2});
global_mid_temp = str2num(A{3});
global_max_temp = str2num(A{4});

% Temperature range
T = linspace(global_min_temp,global_max_temp,500);

% To skip the unwanted lines
for i = 1:3
    fgetl(f1);
end

for i = 1:53
    
    P = fgetl(f1);
    
    Element = strtok(P);
    
    % line-1,2 are required for all high temp.coeffs.
    line_1 = fgetl(f1);
    line_2 = fgetl(f1);
    
    % To find the position of "E"
    a = strfind(line_1,'E');
    b = strfind(line_2,'E');
    
    % High temperature Coefficients
    a1 = str2num(line_1(1:a(1)+3));
    a2 = str2num(line_1(a(1)+3:a(2)+3));
    a3 = str2num(line_1(a(2)+4:a(3)+3));
    a4 = str2num(line_1(a(3)+4:a(4)+3));
    a5 = str2num(line_1(a(4)+4:a(5)+3));
    a6 = str2num(line_1(1:b(1)+3));
    a7 = str2num(line_2(b(1)+4:b(2)+3));
    
    % line-3 is also required for low temp. coeffs.
    line_3 = fgetl(f1);
    c = strfind(line_3,'E');
    
    % Low temperature Coefficients
    a8 = str2num(line_2(b(2)+4:b(3)+3));
    a9 = str2num(line_2(b(3)+4:b(4)+3));
    a10 = str2num(line_2(b(4)+4:b(5)+3));
    a11 = str2num(line_3(1:c(1)+3));
    a12 = str2num(line_3(c(1)+4:c(2)+3));
    a13 = str2num(line_3(c(2)+4:c(3)+3));
    a14 = str2num(line_3(c(3)+4:c(4)+3));
    
    % To calculate Specific heat
    c_p = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,R,T,global_mid_temp);
    
    % To calculate Enthalpy
    H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,R,T,global_mid_temp);
    
    % To calculate Entropy
    S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,R,T,global_mid_temp);
    
    % Plots
    plotting(c_p,H,S,T,Element);
    
    % molecular weight
    mol_wt(i) = molecular_weight(Element);
    disp(mol_wt(i));
    Molecular_Weight = strcat(mol_wt(i),Element);
    
    % write to file [MW]
    Molecular_Weight = mol_wt(i);
    writefile(Element,Molecular_Weight);

end

fclose(f1)

Explanation –

First of all, for every different task, we need to write the function program which we can call later in the main program.
For calculation of specific heat, enthalpy and entropy values we write separate function program with names ‘specific_heat’, ‘enthalpy’ and ‘entropy’ respectively. We use if-else loop in function program of these properties with variation about the global mid temperature (1000).
Then for plotting these thermodynamic properties one by one for all the elements, we create a function named ‘plotting’. Since we want plots of every element in different folder so we use ‘pwd’, ’mkdir’ and ‘cd’ command. ‘pwd’ command is to identify the current dictionary. ‘mkdir’ command is used to make new dictionary. ‘cd’ command is used to change the working dictionary. Then we plot the Cp, H and S plots one by one.
After that, we write the code for the calculation of molecular weights of all elements using multiple for loops and if-else loop. ‘strcmp’ command is used to compare the name of the element with the string array. ‘str2num’ command is used to convert character array or string to numeric array.
As we want a file with the molecular weights of all elements, we create a function ‘writefile’ to create and write the file with molecular weights.
Then we write the main program to calculate everything at once automatically.
First we address the value of the Universal Gas Constant R = 8.314
Then we use the ‘fopen’ command to open ‘THERMO.dat’ file in the read mode specified by 'r' keyword in the ‘fopen’ command.
The file is read line by line using the "fgetl" command. The first line represents the title of the file as “Thermo”
Next step is to obtain the temperature range from the ‘Thermo.dat’ file. The next line has the global maximum, global minimum and global mid values which is read by using 'fgetl' command and split into three values by using the 'strsplit' command and then we have to use ‘str2num’ for converting string cell to numbers.
The temperature array is defined using 'linspace' command as 500 values between the global minimum and global maximum values.
The next three lines are comment in the file so they are ignored by using 'for loop' to iterate for 3 values.
Then 'for' loop is created to loop through the 53 species data given in the file and calculated the required values.
Inside the loop the 'fgetl' command is used to read the next line of the file which represent the name of the species and this variable is given the name ’P’. Function ‘strtok’ is used to define the Elements.
Then we use function ‘strfind’ to find loop inside the loop and use it to derive the value of constants from a1 to a14.
Then we write code to call all the other 6 function to calculate the value of different thermodynamic properties (Cp, H and S), plot these properties with respect to temperature, calculate the molecular weight of each compound and write a file with molecular weight of all the compounds inside the file ‘Molecular_Weight’.

Result –

Value of co-efficients –

File with molecular weights –

Folder with all the plots –

Plots for O2 –

Cp vs Temperature (T) curve -

H vs Temperature (T) curve –

S vs Temperature (T) curve –

Plots for N2 –

Cp vs Temperature (T) curve –

H vs Temperature (T) curve –

S vs Temperature (T) curve –

Plots for CO2 –

Cp vs Temperature (T) curve –

H vs Temperature (T) curve –

S vs Temperature (T) curve –

Reference –

Parsing NASA thermodynamic data using MATLAB

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.