Parsing NASA thermodynamic data using MATLAB

AIM: Write a code in MATLAB to parse NASA thermodynamic data file for various elements.

OBJECTIVES:

1. Write a code that extracts the 14 co-efficient and calculates the enthalpy, entropy and specific heats for all the species in the data file.
2. Calculate the molecular weight of each species and display it in the 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. All these should be generated automatically.
5. Generate Cp, H and S plots.

GOVERNING EQUATION:

Where,

Cp is Specific heat

H is Enthalpy

S is Entropy

R is the universal gas constant = 8.314 J/(mol-K)

T is the Temperature Range

a1, a2, a3, a4, a5, a6, and a7 are the numerical coefficients provided in NASA thermodynamic files.

INPUTS: NASA thermodynamic data file. Link

STEP BY STEP EXPLAINATION:

Before writing the main code for the program we have to write several function codes which are very essential for the working of the main program.

First, we write the code for specific heat function. Here we define the function then we compare the temperature with the global temperature. If temperature is less than global temperature then minimum coefficients are taken into consideration from datafile else maximum coefficients are considered.

Similarly, we write codes for Enthalpy and Entropy.

After that we write code for molecular weight function. Here some elements are defined with their atomic weight then species name for main program are taken as input in the form of string. This string is converted into individual array & compared if species name match with their atomic element. Then we calculate molecular weight using a logic, the logic is to use coefficient number of each species and their atomic weight to calculate molecular weight. This is done by applying for loop. The molecular weight is then printed as output.

Then we write a code which basically just open the molecular weight file, writes the molecular weight of species which is in current loop, save it and close it.

Then we write a code for plotting specific heat, enthalpy and entropy of each species and save it the respective species folder which we create using mkdir command. For each species of element, a separate folder is created with the species name and their plots are saved in that folder only. To achieve this, we make use of mkdir, cd and pwd command. The mkdir command makes a new directory, cd command changes path to that directory and pwd command identifies current directory.

Now, we write the main program code. Initially we open thermo.dat file, fetch the line from the file to use it in the program. Fgetl() command is used to extract data line form the data file and feed it in the program. By using this command, we feed data line by line to the program. The data file contains global temperate range of all gas species which we separate by using strsplit command to split them into global high, medium and low temperature range.

In data file the next three lines are just comments which we skip by using simple for loop. The remaining data is the valuable data containing temperature coefficient for 53 gas species, we extract this data using for loop.

The first line contains species name with the local temperature range. We use strsplit command to split the data specifying that the first value is a string (species name) and other are temperature range which we then assign in low, med and high category according to their position in line. Here the first seven co-efficient are of high temp and other seven are of low temp, we extract them line by line using fgetl() command. After that these co-efficient are used to calculate specific heat, enthalpy and entropy and the plots are made. All this is done with the support of functions that we previous wrote.

CODE:

For Specific heat,

% Specific Heat Function
function CP = SpecificHeat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp)
% Comparing with global temperature
% If less than global it will take minimum coefficients
if Temp > global_medium_temp
CP =(a1+a2.*Temp+a3.*Temp.^2+a4.*Temp.^3+a5.*Temp.^4).*R;
% If temperature is more it will take maximum coefficients
else
CP =(a8+a9.*Temp+a10.*Temp.^2+a11.*Temp.^3+a12.*Temp.^4).*R;
end

For Enthalpy,

% Enthalpy Function
function H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,Temp,R,global_medium_temp)
% Comparing with global temperature
% If less than global it will take minimum coefficients
if Temp > global_medium_temp
H =(a1.*Temp+(a2.*Temp.^2)./2+(a3.*Temp.^3)./3+(a4.*Temp.^4)./4+(a5.*Temp.^5)./5+a6).*R;
else
% If temperature is more it will take maximum coefficients
H =(a8.*Temp+(a9.*Temp.^2)./2+(a10.*Temp.^3)./3+(a11.*Temp.^4)./4+(a12.*Temp.^5)./5+a13).*R;
end

For Entropy,

% Function for entropy
function S = Entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,Temp,R,global_medium_temp)
% Comparing with global temperature
% If less than global it will take minimum coefficients
if Temp > global_medium_temp
S =(a1.*log(Temp)+a2.*Temp+a3.*(Temp.^2)./2+a4.*(Temp.^3)./3+a5.*(Temp.^4)./4+a7).*R;
else
% If temperature is more it will take maximum coefficients
S =(a8.*log(Temp)+a9.*Temp+a10.*(Temp.^2)./2+a11.*(Temp.^3)./3+a12.*(Temp.^4)./4+a14).*R;
end

For Molecular Weight,

% 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

To read and write in text file,

% Function to Read & Write Molecular text file
function Read_write (Species,Molecularweight)
Molecular_list = fopen('Molecular weight.txt','a+');
fprintf( Molecular_list,'Molecular weight of %s is %d \n',Species,Molecularweight);
fclose(Molecular_list);
end

For Plotting,

% For plotting for each species
function plotting(CP,H,S,Temp,Species)
cur_dir=pwd;
mkdir(Species);
cd(Species);
 
% Plot 1, Specific heat (CP) vs Temperature Range (Temp)
figure(1)
plot(Temp,CP,'linewidth',2,'color','r');
xlabel('Temperature [K]');
ylabel('Specific Heat [kJ]');
title(sprintf('Specific Heat vs Temperature Range of %s',Species));
grid on
saveas(figure(1),'Specific Heat.jpg');
 
% Plot 2, Enthalpy (H) vs Temperature Range (Temp)
figure(2)
plot(Temp,H,'linewidth',2,'color','g');
xlabel('Temperature [K]');
ylabel('Enthalpy in [kJ/(mol)]');
title(sprintf('Enthalpy vs Temperature Range of %s',Species));
grid on
saveas(figure(2),'Ethalpy.jpg');
 
% Plot 3, Entropy (S) vs Temperature Range (Temp)
figure(3)
plot(Temp,S,'linewidth',2,'color','b');
xlabel('Temperature [K]');
ylabel('Entropy in [kJ/(mol-K)]');
title(sprintf('Entropy vs Temperature Range of %s',Species));
grid on
saveas(figure(3),'Entropy.jpg');
cd(cur_dir);
 
% Molecular weight
M = Molecular_weight(Species)
f2 = fopen('molecular mass','w');
fprintf('%s n %s n %f n','species','Molecular mass',M)
fclose(f2)
end

Main Code,

clear all
close all
clc
 
% Setting default folder on every startup of Matlab
cd 'C:\Users\sai\Documents\Mat codes\Project'
...
pwd
 
% Opening & Reading NASA Thermodynamic Data file
f1 = fopen ('THERMO.dat','r'); 
gt = fgetl(f1);
% Fetching Header Information
global_temp = fgetl (f1);
temp_range = ('300.000 1000.000 5000.000'); 
% Spliting
A = strsplit(temp_range,' ');
% Converting string to numbers using str2num command
global_low_temp = str2num (A{1}); % Global Low Temp
global_medium_temp = str2num (A{2}); % Global Medium Temp
global_high_temp = str2num (A{3}); % Global High Temp
% skipping the next three comments in the file
for i=3:5
comments = fgetl(f1);
end
for i = 1:53
O_RL = fgetl(f1); % RL = Read line
species = strtok(O_RL);
local_temp = strfind (O_RL,'G');
% reading local temperatue range
Avg_temp = O_RL(local_temp+1:end-1);
Avg_temp = str2num(Avg_temp);
% Range of temperature & Gas Constant
Temp = linspace(global_low_temp,global_high_temp,1000); % Temperature Range
R = 8.314; % Universal Gas Constant J/(mol-K)
% Extracting higher and lower temperature coefficients
line_1 = fgetl(f1);
line_2 = fgetl(f1);
line_3 = fgetl(f1);
A = strfind(line_1,'E');
B = strfind(line_2,'E');
C = strfind(line_3,'E');
% Extracting HIGH-temperature coefficients - FIRST 7 coefficients
a1 = str2num(line_1(1:A(1)+3));
a2 = str2num(line_1(A(1)+4: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_2(1:B(1)+3));
a7 = str2num(line_2(B(1)+4:B(2)+3));
% Extracting LOW-temperature coefficients - SECOND 7 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));
% Specific heat calculation
CP = SpecificHeat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp);
% Enthalpy calculation
H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,Temp,R,global_medium_temp);
% Entropy calculation
S = Entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,Temp,R,global_medium_temp);
% Calculating molecular weight of species
Molecularweight(i) = Molecular_weight(species);
% Read & Write Molecular list
Read_write(species,Molecularweight)
% Creating individual folder for each species and plotting
plotting(CP,H,S,Temp,species)
end

OUTPUT/RESULT/CONCLUSION:

Plots for O₂

Specific Heat

Enthalpy

Entropy

Plots for N₂

Specific Heat

Enthalpy

Entropy

Plots for CO₂

Specific Heat

Enthalpy

Entropy

Command Window,

ERRORS:

1. Silly typo-errors were encountered. Function names are case sensitive in MATLAB so the problem was solved by correctly writing the name of the function.

2. Folder for each gas species were not generating only graphs were generated which were not getting saved. This error was overcome by using pwd, mkdir, cd command in plotting program. As nothing was generated there is no screenshot to show.

REFERENCES, IF ANY:

https://in.mathworks.com/help/matlab/ref/strfind.html

https://in.mathworks.com/help/matlab/ref/strtok.html

https://in.mathworks.com/help/matlab/ref/str2num.html

https://in.mathworks.com/help/matlab/ref/mkdir.html

https://in.mathworks.com/help/matlab/ref/pwd.html

https://in.mathworks.com/help/matlab/ref/cd.html