Parsing NASA Thermodynamic data
NASA Glenn Computer program CEA(Chemical Equilibrium with Applications) is a project which documents thermodynamic data of over 2000 species.This program reported polynomial functions to compute certain thermodynamic properties like Specific Heat, Enathalpy, Entropy, etc. The NASA Glenn Program also compiled a library of data consisting of coefficients to these polynomials specific to species over a wide temperature range. A sample of such data file can be accessed with the following link:
http://combustion.berkeley.edu/gri_mech/version30/files30/thermo30.dat
The aim of this project is to extract these coefficients data from the original dataset and compute thermdynamic parameters Specific Heat(Cp), Enthalpy(H) and Entropy(S) for each of the species. The final putput tabulates molecular weight and coefficients of each species in a .txt file. Also plots for Cp, H and S are automatically created and stored in specific folders named after the species themselves.
As per NASA Glenn Computer program, set of polynomials prescribed for Specific Heat (Cp), Enthalpy (H) and Entropy (E) are as follows:
`(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^-1`
`S/R = a1.lnT + a2.T + a3.T^2/2 + a4.T^3/3 + a5.T^4/4 + a7`
where a1 to a7 are coefficients of polynomials, T the tmeperature in Kelvin scale and R the Gas constant 8.314J/K.
The original dataset contains a set of 14 coefficients out of which first 7 coefficients are high temperature coefficients and rest 7 are low temperature coefficients. In the following program the coefficients are used accordingly in respective temperature ranges to calculate Cp, H and S.
Below is matlab program to extract data from this file and create plots for Specific Heat, Enthalpy and Entropy of the prescribed species and corresponding molecular weight.
The code consists of a main program and five functions. Following five functions have been called to carry out set of operations:
- compute_no_lines: This function computes total number of lines in the original dataset and ouptuts the results to a variable total_lines.
- coff_temp: Extracts coefficients of a specie and respective low, mid and high temperature values from the dataset and stores in a variable data, a column matrix.
- mol_mass: Computes molecular weight of a specie and stores in a variable weight.
- plots_pt: Computes Cp, H and S for a specie over low and high temperature ranges. A folder is created with specie name and plots for Cp, H and S are stored in the folder.
- fnl: Finally this function calls functions 2, 3 and 4, extracts data from them and stores the data in specific variables. This is iterated for all the species in the given dataset in which a tabulation of species-wise molecular weight and coefficients are written to data.txt file.
Below is the main program:
close all
clear all
clc
%Open the data set, read the file and the assign the file to r1
r1 = fopen('thermo30.txt','r');
%Compute number of lines in the file
num_of_lines = compute_no_lines();
%Index each line and assign it to column matrix "l"
for i = 1:num_of_lines
lne = fgetl(r1);
l{i,1} = lne;
end
%Store the location of current directory to path_11
path_11 = pwd
%Create new folder by name "Plots" to store extracted data
mkdir Plots
%Change the directory from current folder to newly created folder "Plots"
cd Plots
%Assign initial value to variable "s" to be used in other functions
s = 7;
%Call function "fnl" to carry out extraction of data and obtain plots
%Compute Molecular Weight
%Extract Coefficients from the original dataset
%Extract low, mid and high temperatures
%Compute Specific Heat, Enthalpy and Entropy over low and high temperatures
%Plot Specific Heat, Enthalpy and Entropy against temperature
%Finally save the individual species molecular wt and coefficients data in a .txt file
plotting = fnl(l,s);
%Return to the earlier directory
Plots = cd(path_11);1. Code for compute_no_lines function
function total_lines = compute_no_lines()
%Open the given data file and read
sd = fopen('thermo30.txt','r');
no_lines = 1;
%fgetl command to extract first line and store it in "s1"
s1 = fgetl(sd);
%Assign data of las line to "s2"
s2 = 'END';
%Compare "s1" and "s2" and iterate untill they are equal
while strcmp(s1,s2) ~= 1
no_lines = no_lines + 1;
%Using fgetl to eaxtract the data line by line
s1 = fgetl(sd);
end
total_lines = no_lines;
end
2. Code for coff_temp function
function data = coff_temp(l,q)
%arguments "l" and "q" obtained fro the function "fnl"
r = 1;
for k = q:(q + 2)
%Extract the coefficients of a specie
c = char(l(k,1));
j = findstr(c,'E');
aa = [];
if k ~= (q + 2)
for m = 1:5
cfft = c((j(m) - 11):(j(m) + 3));
format long
cfft_num(r,1) = str2num(cfft);
r = r + 1;
end
else
for n = 1:4
cffl = c((j(n) - 1):(j(n) + 3));
format long
cfft_num(r,1) = str2num(cffl);
r = r + 1;
end
end
end
coff_a = [];
k = 0;
%Assign the extracted data to a column matrix "coefficients"
coefficients = cfft_num;
%Extract low, mid and high temperature data of corresponding specie
k = q - 1;
lke = char(l(k,1));
low_temp = str2num(lke(49:55));
high_temp = str2num(lke(58:65));
mid_temp = str2num(lke(68:75));
temperature = [low_temp;mid_temp;high_temp];
%Assign the extracted data of coefficients and temperatures to column matrix "data"
format long
data = [coefficients;temperature];
end
3. Code for mol_mass function
function weight = mol_mass(fd)
%argument "fd" obtained from "fnl" funtion
m = 0;
%Atomic mass of C, O, H, N, S and Ar in atomic mass units
mass_C = 12.0107;
mass_O = 15.999;
mass_H = 1.00794;
mass_N = 14.0067;
mass_S = 32.065;
mass_Ar = 39.948;
for i = 1:length(fd)
if fd(i) == 'C'
m = m + mass_C;
elseif fd(i) == 'O'
m = m + mass_O;
elseif fd(i) == 'H'
m = m + mass_H;
elseif fd(i) == 'N'
m = m + mass_N;
elseif fd(i) == 'A'
m = m + mass_Ar;
elseif fd(i) == 'R'
m = m + 0;
elseif fd(i) == '('
m = m + 0;
elseif fd(i) == 'S'
m = m + mass_S;
elseif fd(i) == ')'
m = m + 0;
else
mu = str2num(fd(i));
if fd(i-1) == 'C'
m = m + (mu - 1)*mass_C;
elseif fd(i-1) == 'O'
m = m + (mu - 1)*mass_O;
elseif fd(i-1) == 'H'
m = m + (mu - 1)*mass_H;
else
m = m + (mu -1)*mass_N;
end
end
%Molecular Weight of the species assigned to "weight"
weight = m;
end
4. Code for plots_pt function
function plot_data = plots_pt(data,l,s)
%argumets "data", "l" and "s" obtained from "fnl" function
sc = data;
r = 8.314;
%Create low and high temperatures range to compute Cp, H and S
tl = linspace(sc(15),sc(16),100);
th = linspace(sc(16),sc(17),100);
for i = 1:length(tl)
%Using first eight coefficients for hig temperature range and rest for low temperature range
cp1(i) = (sc(8) + sc(9)*tl(i) + sc(10)*tl(i)^2 + sc(11)*tl(i)^3 + sc(12)*tl(i)^4)*r;
cp2(i) = (sc(1) + sc(2)*th(i) + sc(3)*th(i)^2 + sc(4)*th(i)^3 + sc(5)*th(i)^4)*r;
%Using first eight coefficients for hig temperature range and rest for low temperature range
h1(i) = (sc(8)*tl(i)+(sc(9)*tl(i)^2)/2 + (sc(10)*tl(i)^3)/3 +(sc(11)*tl(i)^4)/4 + (sc(12)*tl(i)^5)/5 + sc(13))*r;
h2(i) = (sc(1)*th(i) + (sc(2)*th(i)^2)/2 + (sc(3)*th(i)^3)/3 +(sc(4)*th(i)^4)/4 + (sc(5)*th(i)^5)/5 + sc(6))*r;
%Using first eight coefficients for hig temperature range and rest for low temperature range
s1(i) = (sc(8)*log(tl(i))+sc(9)*tl(i) + (sc(10)*tl(i)^2)/2 + (sc(11)*tl(i)^3)/3 + (sc(12)*tl(i)^4)/4 + sc(14))*r;
s2(i) = (sc(1)*log(th(i)) + sc(2)*th(i) + (sc(3)*th(i)^2)/2 + (sc(4)*th(i)^3)/3 + (sc(5)*th(i)^4)/4 + sc(7))*r;
end
%Concatenate temperatures, Cp, H and S values to single matrix
temps = [tl,th];
cp = [cp1,cp2];
h = [h1,h2];
st = [s1,s2];
%Create a folder by the name of species to store plot data
nm = char(l((s-1),1));
name = nm(1:7);
mkdir (name)
path_1 = pwd;
%Change directory to species' folder created above
path_2 = pwd;
cd (name)
%Plot Cp and save the file in .png format
figure(1)
plot(temps,cp,'linewidth',2,'color','r')
grid on
xlabel('Temperature(K)')
ylabel('Specific Heat(J/K)')
title(name)
saveas(figure(1),'Specific Heat.png')
%Plot H and save the file in .png format
figure(2)
plot(temps,h,'linewidth',2,'color','b')
grid on
xlabel('Temperature(K)')
ylabel('Enthalpy(J/Kg)')
title(name)
saveas(figure(2),'Enthalpy.png')
%Plot S and save the file in .png format
figure(3)
plot(temps,st,'linewidth',2,'color','g')
grid on
xlabel('Temperature(K)')
ylabel('Entropy(J/K)')
title(name)
saveas(figure(3),'Entropy.png')
%After saving three plot files return to "Plots" directory
Plots = cd(path_2);
plot_data = [temps;cp;h;st];
end
5. Code for fnl function
function data_ex = fnl(l,s)
%arguments "l" and "s" are picked up from main function
%Create a new .txt file to tabulate extracted data
cf = fopen('data.txt','w')
fprintf(cf,'%10s %20s %40s\n\n','Species','Molecular Weight','Coefficients');
%Iteration for 53 species in the original dataset
for e = 1:53
t = s - 1;
%To extract the name of the species
n_m = char(l(t,1));
y = findstr(n_m,' ');
n_n = n_m(1:(y(1)-1));
%Name of species stored in n_n
molecule = n_n;
%Calling the function "mol_mass" to compute molecular weight
molecular_weight = mol_mass(n_n);
%Calling the function "coff_temp" to store coefficients and temperatures in "data"
format long
data = coff_temp(l,s);
%Function "plots_pt" to compute and plot Cp, H and S against %temperatures
plots = plots_pt(data,l,s);
%Iterating for successive species in the original dataset
s = s + 4;
%Converting species name and molecular weight to stirngs to print in the data.txt file
name(e,1) = string(n_n);
mol_weight(e,1) = string(molecular_weight);
%Iteration to extract coefficients from the above created column matrix "data"
for g = 1:14
format long
%Printing Species name, molecular weight and 14 coefficients in order
fprintf(cf,'%10s %20s %40.25f',name(e,1),mol_weight(e,1),data(g));
fprintf(cf,'\n');
end
fprintf(cf,'\n');
end
fprintf(cf,'\n');
fclose(cf);
data_ex = [name,mol_weight];
end
Below are the Cp, H and S plots for O2
Below is the Cp, H and S plots for CO2
Below is the Cp, H and S plots for N2