NASA Thermodynamic Data Parsing in MATLAB

NASA THERMODYNAMIC DATA PARSING

Introduction

NASA came up with ‘Polynomials’ to present thermodynamic properties such as Specific heat, Enthalpy, Entropy as a function of Temperature. Also NASA documented 14 co-efficients for each species to evaluate these polynomials

Inputs

Universal gas constant R = 8.314 J/mol*K

Temperature range = 300K to 5000K

z = pwd;

• ‘pwd’ is used to know current directory.

Open File

Data to be used is present in ‘THERM.dat’ file. To open ‘THERM.dat’ file ‘fopen ‘ command is used, to read data from file ‘r’ is used as permission.

Header Information

Header information mainly contains ‘global minimum temperature, ‘global mid temperature’, ‘global maximum temperature’. It also contains several comments. Steps to extract global minimum, mid and maximum temperature are mentioned below,

• Temperatures are mentioned in second line of ‘THERM.dat’ file. Reach to second line using ‘fgetl’ command and store it in variable ‘tline’.
• As whole line got saved as a single string, ‘strsplit’ command is used to split string and defined ‘ ‘ (space) as delimiter.
• Converted individual temperature string to number by using ‘str2num’ command.
• Next three lines are comments. To skip that lines ‘fgetl’ command is used in ‘for’ loop.

Species Name

After header information, next line contains name of species. Steps to extract species names are mentioned below,

• Reach to line by using ‘fgetl’ command and store it in variable ‘tline’.
• As whole line got saved as a single string, ‘strsplit’ command is used to split string and defined ‘ ‘ (space) as delimiter.
• First string contains species name, to extract species name ‘char(N{1})’ is used and assigned variable ‘name’ to it. Use of ‘strsplit’ command in this condition creates cell array because of that to extract species name ‘{ }’ is used.
• To make folder of every species in parent folder ‘NASA THERMO DATA’ ‘mkdir’ command is used.
• ‘for’ loop is used to make one species’ folder after each loop iteration.

Molecular Weight

Line containing species also contains chemical formula of species (symbols (all capitals) and numbers) after ‘Reference Data Code’. Some compounds in data are made with maximum 4 elements. First element and number of atoms of that element are mentioned in 25^th and 29^th column respectively. Second element and number of atoms of that element are mentioned in 30^th and 34^th column. Third element and number of atoms of that element are mentioned in 35^th and 39^th column respectively. Fourth element and number of atoms of that element are mentioned in 40^th column and 44^th column respectively. To assign molecular weight to particular element comparison of string mentioned as element is done using ‘==’ operator. For comparison ‘if else’ statement is used. To assign molecular weight to element1 below mentioned code is used. This process is repeated 3 more times to assign molecular weight to element2, element3 and element4.

element1 = tline(25);
if element1 == \'H\'
   element1 = 1.00794;
elseif element1 == \'C\'
       element1 = 12.011
elseif element1 == \'N\'
       element1 = 14.00674;
elseif element1 == \'O\'       
       element1 = 15.9994;
elseif element1 == \'A\'
       element1 = 39.948;
end

Code for finding ‘Molecular Weight’ of compound,

if element2 == \' \'
   molecular_weight = (element1*nu_of_atom_element1)
elseif element3 == \' \'
   molecular_weight = (element1*nu_of_atom_element1)+(element2*nu_of_atom_element2)
elseif element4 == \' \'
   molecular_weight = (element1*nu_of_atom_element1)+(element2*nu_of_atom_element2)+(element3*nu_of_atom_element3)
else
   molecular_weight =(element1*nu_of_atom_element1)+(element2*nu_of_atom_element2)+(element3*nu_of_atom_element3)+ (element4*nu_of_atom_element4)

• From above code molecular weight of species can be found according to number of elements it contains.

Steps for Saving Molecular Weight

• Change directory to name. To change directory ‘cd’ command is used.
• Make a text file by using ‘fopen’ command with permission ‘w’ to store molecular weight of each species in it.
• Change directory to current directory.
• By using all code in ‘for’ loop molecular weight of one species is saved in name folder after each loop iteration.

Extracting Co-efficients

To extract coefficients string ‘E’ is found using ‘findstr’ command. Relative to position of string ‘E’ in each co-efficient all co-efficients present in one line can be extracted. To convert co-efficients to number from strings ‘str2num’ command is used. To extract all 14 co-efficients of one species this process is repeated for next lines as co-efficients are mentioned in 3 lines. To extract co-efficients for all species ‘for’ loop is used.

Specific heat, Enthalpy, Entropy

To calculate Specific heat, Enthalpy and Entropy individual functions are made as per the formulas. Plots are made using ‘plot’ command. Steps to save plots are mentioned below,

• Change directory to name using ‘cd’ command.
• Save plot in name folder in ‘jpg’ format using ‘saveas’ command.
• Change directory to current directory.

Main Program

clear all
close all
clc

%Inputs
%Universal gas constant in J/(mol*k)
R = 8.314;
T = linspace(300,5000,200);
number_of_species = 53;
z = pwd;

f1 = fopen(\'THERMO.dat\',\'r\');
tline = fgetl(f1);
tline = fgetl(f1);
A = strsplit(tline,\' \');
global_low_temp = str2num(A{2});
global_mid_temp = str2num(A{3});
global_high_temp = str2num(A{4});

%Header Information
for i = 1:3;
    tline=fgetl(f1);
end
%Species Name
for j=1:number_of_species
tline = fgetl(f1);
N = strsplit(tline,\' \');
name = char(N{1});
mkdir(\'NASA THERMO DATA\',name);
%Molecular Weight of Species
element1 = tline(25);
if element1 == \'H\'
   element1 = 1.00794;
elseif element1 == \'C\'
       element1 = 12.011;
elseif element1 == \'N\'
       element1 = 14.00674;
elseif element1 == \'O\'       
       element1 = 15.9994;
elseif element1 == \'A\'
       element1 = 39.948;
end
nu_of_atom_element1 = tline(29);
nu_of_atom_element1 = str2num(nu_of_atom_element1);
element2 = tline(30);
if element2 == \'H\'
   element2 = 1.0074;
elseif element2 == \'C\'
       element2 = 12.011;
elseif element2 == \'N\'
       element2 = 14.00674;
elseif element2 == \'O\'       
       element2 = 15.9994;
elseif element2 == \'A\'
       element2 = 39.948;
end
nu_of_atom_element2 = tline(34);
nu_of_atom_element2 = str2num(nu_of_atom_element2);
element3 = tline(35);
if element3 == \'H\'
   element3 = 1.0074;
elseif element3 == \'C\'
       element3 = 12.011;
elseif element3 == \'N\'
       element3 = 14.00674;
elseif element3 == \'O\'       
       element3 = 15.9994;
elseif element3 == \'A\'
       element3 = 39.948;
end
nu_of_atom_element3 = tline(39);
nu_of_atom_element3 = str2num(nu_of_atom_element3);
element4 = tline(40);
if element4 == \'H\'
   element4 = 1.0074;
elseif element4 == \'C\'
       element4 = 12.011;
elseif element4 == \'N\'
       element4 = 14.00674;
elseif element4 == \'O\'       
       element4 = 15.9994;
elseif element4 == \'A\'
       element4 = 39.948;
end
nu_of_atom_element4 = tline(44);
nu_of_atom_element4 = str2num(nu_of_atom_element4);
if element2 == \' \'
   molecular_weight = (element1*nu_of_atom_element1);
elseif element3 == \' \'
   molecular_weight = (element1*nu_of_atom_element1)+(element2*nu_of_atom_element2);
elseif element4 == \' \'
   molecular_weight = (element1*nu_of_atom_element1)+(element2*nu_of_atom_element2)+(element3*nu_of_atom_element3);
else
   molecular_weight =(element1*nu_of_atom_element1)+(element2*nu_of_atom_element2)+(element3*nu_of_atom_element3)+ (element4*nu_of_atom_element4);
end
cd(\'NASA THERMO DATA\');
cd(name);
f2 = fopen(\'Molecular Weight.txt\',\'w\');
fprintf(f2,\'Molecular weight of species = %f\', molecular_weight);
cd(z);
%Extracting Higher Temperature Coeffs
tline = fgetl(f1);
C = findstr(tline,\'E\');
a1 = tline(1:C(1)+3);
a1 = str2num(a1);
a2 = tline(C(1)+4:C(2)+3);
a2 = str2num(a2);
a3 = tline(C(2)+4:C(3)+3);
a3 = str2num(a3);
a4 = tline(C(3)+4:C(4)+3);
a4 = str2num(a4);
a5 = tline(C(4)+4:C(5)+3);
a5 = str2num(a5);
tline = fgetl(f1);
C = findstr(tline,\'E\');
a6 = tline(1:C(1)+3);
a6 = str2num(a6);
a7 = tline(C(1)+4:C(2)+3);
a7 = str2num(a7);
%Extracting Lower Temperature Coeffs
a8 = tline(C(2)+4:C(3)+3);
a8 = str2num(a8);
a9 = tline(C(3)+4:C(4)+3);
a9 = str2num(a9);
a10 = tline(C(4)+4:C(5)+3);
a10 = str2num(a10);
tline = fgetl(f1);
a11 = tline(1:C(1)+3);
a11 = str2num(a11);
a12 = tline(C(1)+4:C(2)+3);
a12 = str2num(a12);
a13 = tline(C(2)+4:C(3)+3);
a13 = str2num(a13);
a14 = tline(C(3)+4:C(4)+3);
a14 = str2num(a14);
%Plotting Temperature vs Specific heat
cp = specific_heat(R,T,a1,a2,a3,a4,a5,a8,a9,a10,a11,a12);
h = figure(1);
plot(cp,T,\'linewidth\',2);
xlabel(\'specific heat\');
ylabel(\'Temperature\');
title(name);
cd(\'NASA THERMO DATA\');
cd(name);
saveas(h,\'SPECIFIC HEAT.jpg\');
cd(z);
%Plotting Temperature vs Enthalpy
H = enthalpy(R,T,a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13);
figure(2)
v = figure(2);
plot(H,T,\'linewidth\',2);
xlabel(\'Enthalpy\');
ylabel(\'Temperature\');
title(name);
cd(\'NASA THERMO DATA\');
cd(name);
saveas(v,\'ENTHALPY.jpg\');
cd(z);
%Plotting Temperature vs Entropy
S = entropy(R,T,a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14);
figure(3)
w = figure(3);
plot(S,T,\'linewidth\',2);
xlabel(\'Entropy\');
ylabel(\'Temperature\');
title(name);
cd(\'NASA THERMO DATA\');
cd(name);
saveas(w,\'ENTROPY.jpg\');
cd(z);
end

fclose(f1);

Functions

function [cp] = specific_heat(R,T,a1,a2,a3,a4,a5,a8,a9,a10,a11,a12)

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

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

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

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

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

Images