Parsing NASA thermodynamic data

AIM - To parse the NASA thermodynamic data using MATLAB.

Theory - There is a NASA thermodynamic data having different thermodynamic values w.r.t. different temperatures in polynomial form in '.dat' format also having 14 coefficients. In order to extract the different information we need to parse the data using MATLAB.

The term ‘parsing’ refers to MATLAB reading an m-file (script or function) and translating it into an internal code that it stores in memory. It is this ‘parsed’ file in memory that is actually run when you execute the code.

Governing equations -

Specific heat -

`C_p/R =a_1+a_2*T+a_3*T^2+a_4*T^3+a_5*T^4`

Enthalpy -

`H/(RT) = a_1 +a_2*T/2+a3*T^2/3+a_4*T^3/4+a_5*T^4/5+a_6/T`

Entropy -

`S/R = a_1*lnT+a_2*T+a_3*T^2/2+a_4*T^3/3+a_5*T^4/4 +a_7`

where,

C_P = Specific Heat

H = Enthalpy

S = Entropy

R = Universal Gas Constant (8.314 J/mol.K)

T = Temperature

The First 7 Coefficients are High Temperature Coefficients and the Second 7 Coefficients are Low Temperature Coefficients.

Objectives -

1.Reading "NASA thermodynamic data" file and write a function that extracts the 14 coefficients and calculates the enthalpy, entropy and specific heats for all the species.

2.To calculate molecular weight of every species and displaying in command window.

3.Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species.

4.Save the plots as images with appropriate names and dump them in separate folders for each species with suitable name.

Procedure -

Function for Molecular weight calculation for different species-

% Function for Molecular Weight
function Molecularweight = Molecular_weight(Species)
atomic_name = ['H','C','O','N','A']; 
fprintf('Molecular Weight of ');
fprintf(Species);
fprintf(' : ');
Atomic_Weight = [1.008 12.011 15.999 14.007 39.948];
Molecularweight = 0; % Initiating 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

Main Program -

clear all
close all
clc

f1 = fopen('THERMO.dat','r');
fgetl(f1);

% Defining temperature values
templ = str2num(fgetl(f1));
global_low_temp = templ(1);
global_mid_temp = templ(2);
global_high_temp = templ(3);

% For skipping unwanted lines
for i = 1:3
fgetl(f1);
end

for i=1:53

line1 = fgetl(f1);
l = strsplit(line1, ' ');
m = length(l);
species = l{1};

% Storing the values of the local temperature values

low_temp = str2double(l(m-3));
high_temp = str2double(l(m-2));
mid_temp = str2double(l(m-1)); % Check this later

% Storing the values of the second line of variables
line2 = fgetl(f1);

% Finding the position of E
a = strfind(line2,'E');

% Reading all the other variables
a1 = str2double(line2(1:(a(1)+3)));
a2 = str2double(line2((a(1)+4):(a(2)+3)));
a3 = str2double(line2((a(2)+4):(a(3)+3)));
a4 = str2double(line2((a(3)+4):(a(4)+3)));
a5 = str2double(line2((a(4)+4):(a(5)+3)));

% Storing the values of the 3rd line
line3 = fgetl(f1);
% Finding the position of E
b = strfind(line3,'E');



% Storing the values for the next set of variables
a6 = str2double(line3(1:b(1)+3));
a7 = str2double(line3((b(1)+4):(b(2)+3)));
a8 = str2double(line3((b(2)+4):(b(3)+3)));
a9 = str2double(line3((b(3)+4):(b(4)+3)));
a10 = str2double(line3((b(4)+4):(b(5)+3)));

% Storing ther values of the 4th line
line4 = fgetl(f1);

% Finding the position of E
c = strfind(line4,'E');

% Storing the values of variables in line 4
a11 = str2double(line4(1:c(1)+3));
a12 = str2double(line4((c(1)+4):(c(2)+3)));
a13 = str2double(line4((c(2)+4):(c(3)+3)));
a14 = str2double(line4((c(3)+4):(c(4)+3)));


% Creating Temperature values for lower and higher temperature
lower_temp = linspace(low_temp,mid_temp,100);
higher_temp = linspace(mid_temp,high_temp,100);

% Value of the universal has constant in j/mol*K
R = 8.314; % Universal Gas Constant


% Calculating the values for cp, H and S using the given formulas
for k=1:100

% lower temperature values for cp,h,s calculations i.e  low to middle temp values
cp_low(k)=(a8+(a9*lower_temp(k))+(a10*(lower_temp(k)^2))+(a11*(lower_temp(k))^3)+(a12*(lower_temp(k))^4))*R;

h_low(k)=(a8+(a9*lower_temp(k))/2+(a10*lower_temp(k)^2)/3+(a11*lower_temp(k)^3)/4+a12*(lower_temp(k)^4)/5+(a13/lower_temp(k)))*R*lower_temp(k);

s_low(k)=((a8*log(lower_temp(k))+(a9*lower_temp(k))+((a10*lower_temp(k)^2)/2)+((a11*lower_temp(k)^3)/3)+((a12*lower_temp(k)^4)/4)+a14)*R);

% higher temperature values for cp,h,s calculations i.e mid to high temp values
cp_high(k)=(a1+(a2*higher_temp(k))+(a3*(higher_temp(k)^2)+(a4*(higher_temp(k))^3)+a5*(higher_temp(k))^4))*R;

h_high(k)=(a1+(a2*higher_temp(k))/2+(a3*higher_temp(k)^2)/3+(a4*higher_temp(k)^3)/4+(a5*(higher_temp(k)^4))/5+(a6/higher_temp(k)))*R*higher_temp(k);

s_high(k)=((a1*log(higher_temp(k))+a2*higher_temp(k)+((a3*higher_temp(k)^2)/2)+((a4*higher_temp(k)^3)/3)+((a5*higher_temp(k)^4)/4)+a7)*R);

end

%Molecular weight calculation
Molecularweight(i) = Molecular_weight(species);
% Making Folders for saving the plots
nasa_graphs = species;
mkdir(nasa_graphs);
current = pwd;
cd(nasa_graphs);

%Plotting the graphs

% For specific heat
figure(1)
plot(lower_temp,cp_low,'linewidth',2,'color','r');
hold on
plot(higher_temp,cp_high,'linewidth',2,'color','r');
hold off
title(['Temperature vs Specific Heat',species]);
xlabel('Temperature');
ylabel('Specific Heat');
legend('low temperature','high temperature');
fig = sprintf('Specific Heat of %s.png',species);
saveas(gca,fig);

% For Enthalpy
figure(2)
plot(lower_temp,h_low,'linewidth',2,'color','c');
hold on
plot(higher_temp,h_high,'linewidth',2,'color','c');
hold off
title(['Temperature vs Enthalpy',species]);
xlabel('Temperature');
ylabel('Enthalpy');
legend('low temperature','high temperature');
fig = sprintf('Enthalpy of %s.png',species);
saveas(gca,fig);

% For Entropy
figure(3)
plot(lower_temp,s_low,'linewidth',2,'color','y');
hold on
plot(higher_temp,s_high,'linewidth',2,'color','y');
hold off
title(['Temperature vs Entropy',species]);
xlabel('Temperature');
ylabel('Entropy');
legend('low temperature','high temperature');
fig = sprintf('Entropy of %s.png',species);
saveas(gca,fig);
cd(current)
end

Code Explanation -

1. First write the molecular weight function to call up it in the main program to display it in command window.

2. Read the NASA thermodynamic data by using fopen(to open) and then 'r' command to read it.

3. Extract the global temperatures and store it from NASA file.

4. To extract the local temperatures from the file.

5. Read all the coefficients of NASA data file from a1 to a14 and convert that string to number using 'str2double'.

6. Make a folder/directory named as named as 'nasa data' in 'C:' to store all the graphs

7. Calculate the specific heat(Cp),Enthalpy(H), Entropy(S) by the thermodynamic relations and molecular weight by calling the function.

8.Plot the datas using plot command and save the fig. in png format using saveas command.

9. Molecular weight will be displayed in command window.

Results and Outputs-

Molecular weights in command window of MATLAB-

Plots of O2 -

Plots of N2 -

Plots of Co2 -

Folder image generated automatically -

Conclusion -

The NASA thermodynamic data file is parsed.

The NASA THERMO.dat file -

https://drive.google.com/open?id=1AikCHDlz9s_fu8whE81qee_GD4AyYpCP