Parsing NASA thermodynamic data using MATLAB

Objective: To write a function that extracts the 14 coefficients from NASA thermodynamic data file and calculates the enthalpy, entropy and specific heats for all the species in the data file

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

The polynomials are given as:

Cp/R = a1 + a2*T + a3*T^2 + a4*T^3 + a5*T^4

H/RT = a1 + a2*T/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 + a6 + a7

MATLAB has the feature to read and write .txt files using fopen command. This is used to read and extract the polynomials from the NASA thermodynamic data file. (a1, a2, a3, ....a14).

The first seven coefficients are for lower temperature range and the last seven coefficients are for higher temperature range.

Procedure:

1. The NASA thermodynamic data file is read using fopen command
2. The unwanted header at the start of the file is skipped by using fgetl command
3. The species name is extracted using fgetl command and storing it in a pointer
4. Since the temperature values contain decimal points, they are extracted by finding '.' in the line.
5. The temperature matrix is defined by considering the temperature range from the extracted values
6. Since the coefficients contain letter 'E', they are found using findstr command and then extracted and stored in the variables
7. Thus, all the coefficients of high and low temperature ranges are obtained for a particular species
8. Using for loop, the coefficients are extracted and stored for all the species in the data file
9. Now, using these extracted coefficients, the properties ssuch as Cp, H & S are calculated using the polynomial formula
10. These properties are plotted and saved in different species folder.
11. The molecular weight is also calculated for each species and displayed in the command window

MATLAB code:

close all
clear all
clc
% Parsing NASA thermodynamic data
f1 = fopen('THERMO.dat','r');
for i=1:5
fgetl(f1);
end


for i=1:53      %for 53 species in data file
n=500;
line1 = fgetl(f1);
T = strsplit(line1,' ');
species = T{1};
a=findstr(line1,'.');
T_lower =str2num(line1(a(1)-3:a(1)+3));
T_middle =str2num(line1(a(2)-4:a(2)+3));
T_higher=str2num(line1(a(3)-4:a(3)+3));
t = linspace(T_lower,T_higher,n);

line2 = fgetl(f1);
n2 = findstr(line2,'E');

%Extracting and saving low temperture coefficients
a1_l = str2num(line2(1:n2(1)+3));
a2_l = str2num(line2(n2(1)+4:n2(2)+3));
a3_l = str2num(line2(n2(2)+4:n2(3)+3));
a4_l = str2num(line2(n2(3)+4:n2(4)+3));
a5_l = str2num(line2(n2(4)+4:n2(5)+3));
line3 = fgetl(f1);
n3 = findstr(line3,'E');
a6_l = str2num(line3(1:n3(1)+3));
a7_l = str2num(line3(n3(1)+4:n3(2)+3));

%Extracting and saving high temperture coefficients
a1_h =str2num(line3(n3(2)+4:n3(3)+3));
a2_h = str2num(line3(n3(3)+4:n3(4)+3));
a3_h = str2num(line3(n3(4)+4:n3(5)+3));
line4 = fgetl(f1);
n4 = findstr(line4,'E');
a4_h = str2num(line4(1:n4(1)+3));
a5_h= str2num(line4(n4(1)+4:n4(2)+3));
a6_h =str2num(line4(n4(2)+4:n4(3)+3));
a7_h = str2num(line4(n4(3)+4:n4(4)+3));

R =  8.314 ;

%Property Calculation
for i=1:n
if  (t (i)<= T_middle)
Cp(i)= R*(a1_l + a2_l*t(i) + a3_l*t(i).^2 + a4_l*t(i).^3 + a5_l*t(i).^4);
H(i)= R*t(i)*(a1_l + a2_l*t(i)/2 + a3_l*(t(i)^2)/3 + a4_l*(t(i)^3)/4 + a5_l*(t(i)^4)/5 + a6_l/t(i));
S(i) = R*(a1_l*log(t(i)) + a2_l*t(i) + a3_l*(t(i)^2)/2 + a4_l*(t(i)^3)/3 + a5_l*(t(i)^4)/4 + a7_l);
else 
Cp(i)= R*(a1_h + a2_h*t(i) + a3_h*t(i)^2 + a4_h*t(i)^3 + a5_h*t(i)^4);
H(i) = R*t(i)*(a1_h + a2_h*t(i)/2 + a3_h*(t(i)^2)/3 + a4_h*(t(i)^3)/4 + a5_h*(t(i)^4)/5 + a6_l/t(i));
S(i) = R*(a1_h*log(t(i)) + a2_h*t(i) + a3_h*(t(i)^2)/2 + a4_h*(t(i)^3)/3 + a5_h*(t(i)^4)/4 + a7_h);
endif
end

cur_dir=pwd;
mkdir(species);
cd(species);

figure(1)
plot(t,Cp,'linewidth',3)
title(sprintf('Specific heat capacity vs Temperature for %s',species))
ylabel('Specific Heat capacity');
xlabel('Temperature');
grid on
saveas(figure(1),'Specific heat.png');
pause(0.01)

figure(2)
plot(t,H,'linewidth',3)
title(sprintf('Enthalpy vs Temperature for %s',species))
ylabel('Enthalpy');
xlabel('Temperature');
grid on
saveas(figure(2),'Enthalpy.png');
pause(0.01)

figure(3)
plot(t,S,'linewidth',3)
title(sprintf('Entropy vs Temperature for %s',species))
ylabel('Entropy');
xlabel('Temperature');
grid on
saveas(figure(3),'Entropy.png');
pause(0.01)

cd(cur_dir);

atomic_name = ['H','C','O','N','A'];             % Standard Atomic Name
 % Hydrogen, Carbon, Oxygen, Nitrogen, Argon
 
 fprintf('Molecular Weight of ');
 fprintf(species);
 fprintf(' : ');
 Atomic_Weight = [1.008 12.011 15.999 14.007 39.948];     % Standard Atomic Number
 Molecularweight = 0;                          % Starting Range of Molecule
 % Adding Atomic number to get particular value
 for i = 1:length(species)
 for j = 1:length(atomic_name)
 if strcmpi(species(i),atomic_name(j))
 Molecularweight = Molecularweight +Atomic_Weight(j);
 position = j;
 end
 end
 n = str2num(species(i));
 if n>1
 Molecularweight = Molecularweight + Atomic_Weight(position)*(n-1);
 end
 end
 fprintf('molecular weight %s :',species);
 fprintf('%d',Molecularweight);
 disp(' ');
end

Results:

Command window displaying molecular weight of species:

Plots for Specific Heat, Enthalpy & Entropy for O2:

Plots for Specific Heat, Enthalpy & Entropy for N2:

Plots for Specific Heat, Enthalpy & Entropy for CO2:

Automatic dumping of plots in their respective folders: