MATLAB Project: Parsing NASA thermodynamic data
Aim:
To write a code in MATLAB that extracts the 14 co-efficients and calculates the thermodynamic properties (enthalpy, entropy, molecular weight, and specific heats) of all the species in the data file 'NASA thermodynamic data'.
Theory:
Parsing means to break it down into component parts so that its syntax can be analyzed, categorized, and understood [1]. By doing this we are going to a string into parts then recognizing the parts to convert it into something simpler than a string. File parsing is to resolve the things into component or subpart form and describe the syntactic file. That is, reading in data stream and building a memory model of it.
NASA Thermodynamic File
NASA came up with polynomials that can be used to evaluate the thermodynamic properties such as Enthalpy, Entropy, and Specific heat. They have also documented the co-efficients that are required to evaluate these polynomials for 1000+ species.
NASA thermodynamic file which is 'THERMO.dat' file contains the temperature range (global), 14 coefficients (both low and high temperature coefficients) and the local temperature range of the 53 species.
Figure 1: Sample view of THERMO.dat
Figure 1 shows the sample view of THERMO.dat, in this dat file first line is the heading 'THERMO'. The next line is the global temperature range, where 300.000 is the minimum global temperature while 1000.000 and 5000.000 are medium and maximum global temperature respectively. The next three lines are general information which is not required for the calculation of thermodynamic properties. The next each four lines provides the data of the species which is required for the calculation of thermodynamic properties. Out of these four lines, first line indicates the species and its local temperature range. The next three lines provides the temperature coefficient which is required for the calculation of H, S and Cp. The next four lines are for other species and it is repeated for 53 elements.
The 53 elements in 'THERMO.dat' file are O, O2, H, H2, OH, H2O, HO2, H2O2, C, CH, CH2, CH2(S), CH3, CH4, CO, CO2, HCO, CH2O, CH2OH, CH3O, CH3OH, C2H, C2H2, C2H3, C2H4, C2H5, C2H6, CH2CO, HCCO, HCCOH, H2CN, HCN, HNO, N, NNH, N2O, NH, NH2, NH3, NO, NO2, HCNO, HOCN, HNCO, NCO, CN, HCNN, N2, AR, C3H8, C3H7, CH3CHO, and CH2CHO.
Thermodynamic Properties
Enthalpy, entropy and specific heat of various species at different temperatures are required extensively to analyze the thermodynamic process. For example, if the combustion process is considered, these quantities are required to observe how the combustion is going on and how efficient or effective is. It then allows to seek for the improving the process. Similarly, in other processes like refrigeration, air conditioning, nozzle flow, chemical kinetics, etc, these quantities are widely involved to observe the phenomenon.
Enthalpy
Enthalpy is a property of a thermodynamic system which is equal to the system's internal energy plus the product of its pressure and volume. In a system contained so as to prevent mass transfer, for processes at constant pressure, the heat absorbed or released equals the change in enthalpy.
Entropy
Entropy is the measure of system's thermal energy per unit temperature that is unavailable for doing useful work. Because work is obtained from ordered molecular motion, the amount of entropy is also a measure of the molecular disorder, or randomness of a system.
Specific Heat
Specific heat is the amount of heat per unit mass required to raise the temperature by one degree Celsius.
Molecular Weight
The average mass of a compound compared to 1/12 the mass of carbon 12 and calculated as the sum of the atomic weights of the constituent atoms.
Atomic Mass
Atomic mass is the mass of a single atom of a chemical element. It includes the masses of the 3 subatomic particles that make up an atom: protons, neutrons and electrons. Atomic mass can be expressed in grams.
The empirical equation to find the thermodynamic properties using the coefficients from the THERMO.dat is provided below:
Enthalpy, `H = RT(a_1+a_2T/2+a_3T^2/3+a_4T^3/4+a_5T^4/5+a_6/T)`
Entropy, `S= R(a_1ln(T)+a_2T+a_3T^2/2+a_4T^3/3+a_5T^4/4+a_7)`
Specific Heat, `Cp=R(a_1+a_2T+a_3T^2+a_4T^3+a_5T^4)`
where, Cp is Specific Heat [J/(K kg)], R is Universal Gas Constant (R=8.314 [J/(mol K)]), H is the Enthalpy [J], S is the Entropy [J/K] and T is the temperature [K].
Objectives:
- Understanding the theoretical concept of File Parsing and the commands required for it.
- Writing the program for file parsing.
- Understanding the concepts of thermodynamic properties and the equations required to calculate it.
- Incorporating the functions to calculate Entropy, Enthalpy, Specific Heat and Molecular Weight to the main program and plotting the graphs.
- Save the plots as images with appropriate names and dump them in separate folders for each species with suitable name.
Program:
The programming is done in MATLAB.
Main Code:
% Parsing NASA Thermodynamic data
clear all
close all
clc
% opening the file (NASA Thermodynaic data)
f1 = fopen('THERMO.dat','r'); % for extracting datas
f2 = fopen('Molecular Weights.txt','w'); % for writing molecular weight to new text file
% getting heading
header = fgetl(f1);
% getting global temperature values
global_temp = fgetl(f1);
global_temp = str2num(global_temp);
global_low_temp = global_temp(1);
global_medium_temp = global_temp(2);
global_high_temp = global_temp(3);
% Universal Gas constant [J/(mol-K)]
R = 8.314;
% Skipping next 3 comment lines
for i = 1:3
tline = fgetl(f1);
end
% Extracting coefficients
% Note: 53 species, each species with 14 coefficients (7 high & 7 low temp
% coefficients
for i = 1:53 % 53 species
tline = fgetl(f1);
s = strsplit(tline,' ');
species = (s{1}); % extracting each species name
% Finding local temperature
local_medium_temp = str2num(s{end-1});
local_high_temp = str2num(s{end-2});
local_low_temp = str2num(s{end-3});
% setting time for both lower and higher coeff range // FOR GRAPH
t1 = linspace(local_low_temp,local_medium_temp,1000); % temperature range for low temp coeff
t2 = linspace(local_medium_temp,local_high_temp,1000); % temperature range for high temp coeff
% the coefficients are written in 3 lines, 5 coeffs each in first 2
% lines and 4 coeffs in third line
line1 = fgetl(f1);
line2 = fgetl(f1);
line3 = fgetl(f1);
a = strfind(line1,'E'); % for finding the position of E which is part of every coeff for each line // for line 1
b = strfind(line2,'E'); % for line 2
c = strfind(line3,'E'); % for line 3
% Extracting 7 high temp coeffs
a1 = line1(1:a(1)+3); % coeff 1
a1 = str2num(a1);
a2 = line1(a(1)+4:a(2)+3); % coeff 2
a2 = str2num(a2);
a3 = line1(a(2)+4:a(3)+3); % coeff 3
a3 = str2num(a3);
a4 = line1(a(3)+4:a(4)+3); % coeff 4
a4 = str2num(a4);
a5 = line1(a(4)+4:a(5)+3); % coeff 5
a5 = str2num(a5);
a6 = line2(1:b(1)+3); % coeff 6
a6 = str2num(a6);
a7 = line2(b(1)+4:b(2)+3); % coeff 7
a7 = str2num(a7);
% Extracting 7 low temp coeffs
a8 = line2(b(2)+4:b(3)+3); % coeff 8
a8 = str2num(a8);
a9 = line2(b(3)+4:b(4)+3); % coeff 9
a9 = str2num(a9);
a10 = line2(b(4)+4:b(5)+3); % coeff 10
a10 = str2num(a10);
a11 = line3(1:c(1)+3); % coeff 11
a11 = str2num(a11);
a12 = line3(c(1)+4:c(2)+3); % coeff 12
a12 = str2num(a12);
a13 = line3(c(2)+4:c(3)+3); % coeff 13
a13 = str2num(a13);
a14 = line3(c(3)+4:c(4)+3); % coeff 14
a14 = str2num(a14);
% Entropy Calculation
[S1,S2] = entropy_calc(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,t1,t2,R);
% Enthalpy Calculation
[H1,H2] = enthalpy_calc(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,t1,t2,R);
% Specific Heat Calculation
[Cp1,Cp2] = specificheat_calc(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,t1,t2,R);
% Molecular Weight Calculation
molecular_weight(i) = molecularweight_calc(species);
k = sprintf('Molecular weight of %s is %f u',species,molecular_weight(i));
disp(k);
% Saving Molecular weight in a txt file
fprintf(f2,'Molecular weight of %s is %f u n',species,molecular_weight(i));
% Creating folder
mkdir('C:UsersjibinOneDriveDesktopMATAB programsweek 4Parsing',species);
% Plotting graph
% Temperature vs Entropy
figure(1)
plot(t1,S1,'linewidth',2,'color','b');
hold on;
plot(t2,S2,'linewidth',2,'color','r');
hold off;
xlabel('Temperature [K]');
ylabel('Entropy [J/K]');
title(species);
legend('low temperature coefficient','high temperature coefficient');
grid on;
grid minor;
% Temperature vs Enthalpy
figure(2)
plot(t1,H1,'linewidth',2,'color','b');
hold on;
plot(t2,H2,'linewidth',2,'color','r');
hold off;
xlabel('Temperature [K]');
ylabel('Enthalpy [J]');
title(species);
legend('low temperature coefficient','high temperature coefficient');
grid on;
grid minor;
% Temperature vs Specific Heat
figure(3)
plot(t1,Cp1,'linewidth',2,'color','b');
hold on;
plot(t2,Cp2,'linewidth',2,'color','r');
hold off;
xlabel('Temperature [K]');
ylabel('Specific heat [J/kg/K]');
title(species);
legend('low temperature coefficient','high temperature coefficient');
grid on;
grid minor;
% Saving the figures in the corresponding folder
location = (['C:UsersjibinOneDriveDesktopMATAB programsweek 4Parsing',species]);
cd(location)
saveas(1,'Entropy.jpg');
saveas(2,'Enthalpy.jpg');
saveas(3,'Specific heat.jpg');
cd('C:UsersjibinOneDriveDesktopMATAB programsweek 4');
end
Each step of the program is explained with comments.
Function to calculate Enthalpy:
% Enthalpy Calculation
% H1 and t1 corresponds to the lower temperature range where the
% coefficients a8 to a14 are to be used.
% H2 and t2 corresponds to the higher temperature range where the
% coefficients a1 to a7 are to be used.
function [H1,H2] = enthalpy_calc(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,t1,t2,R)
H2 = R.*t2.*((a1)+(a2.*(t2./2))+(a3.*((t2.^2)./3))+(a4.*((t2.^3)./4))+(a5.*((t2.^4)./5))+(a6./t2));
H1 = R.*t1.*((a8)+(a9.*(t1./2))+(a10.*((t1.^2)./3))+(a11.*((t1.^3)./4))+(a12.*((t1.^4)./5))+(a13./t1));
endFunction to calculate Entropy:
% Entropy calculation function
% S1 and t1 corresponds to the lower temperature range where the
% coefficients a8 to a14 are to be used.
% S2 and t2 corresponds to the higher temperature range where the
% coefficients a1 to a7 are to be used.
function [S1,S2] = entropy_calc(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,t1,t2,R)
S2 = R.*((a1.*log(t2))+(a2.*t2)+(a3.*((t2).^2)/2)+(a4.*((t2).^3)/3)+(a5.*((t2).^4)/4)+(a7));
S1 = R.*((a8.*log(t1))+(a9.*t1)+(a10.*((t1).^2)/2)+(a11.*((t1).^3)/3)+(a12.*((t1).^4)/4)+(a14));
end Function to calculate Specific heat:
% Specific Heat Calculation
% Cp1 and t1 corresponds to the lower temperature range where the
% coefficients a8 to a14 are to be used.
% Cp2 and t2 corresponds to the higher temperature range where the
% coefficients a1 to a7 are to be used.
function [Cp1,Cp2] = specificheat_calc(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,t1,t2,R)
Cp2 = R.*((a1)+(a2.*t2)+(a3.*(t2.^2))+(a4.*(t2.^3))+(a5.*(t2.^4)));
Cp1 = R.*((a8)+(a9.*t1)+(a10.*(t1.^2))+(a11.*(t1.^3))+(a12.*(t1.^4)));
endFunction to calculate Molecular Weight:
% Molecular Weight Calculation
function mol_wt = molecularweight_calc(species)
% Standard Elements present in the data
element = {'H' 'C' 'N' 'O' 'S' 'A'};
% Hydrogen, Carbon, Nitrogen, Oxygen, Argon
% Atomic no. of each elements
atomic_mass = [1.008 12.011 15.999 14.007 32.065 39.948];
mol_wt = 0;
for i = 1:length(species)
for j = 1:length(element)
if strcmp(species(i),element(j))
mol_wt = mol_wt+atomic_mass(j);
pos = j;
end
end
num = str2double(species(i));
if num > 1
mol_wt = mol_wt+(atomic_mass(pos)*(num-1));
end
end
endEach step of the program is explained with comments.
Logic of the program
To get access to 'THERMO.dat', the command fopen is used, and to get(extract) the content of each line, the command fgetl is used. The datas from each line is extracted or separated using the command strsplit and each splitted data are actually in character form. In order to change the form on numerical datas from character, str2num [2] command is used. For calculating the thermodynamic properties the datas required such as 14 coefficients, 3 local temperatures (low, medium and high) are extracted. This is calculated using the formula through different function for calculating each thermodynamic property. In order to show two regions of low and high temperature range in graph, H1, S1 and Cp1 (for lower temperature range) and H2, S2 and Cp2 (for higher temperature range) are calculated. The coefficients are extracted by the logic such that 'E' is found from each line since every coefficient contains E and there are no other Es present. Using mkdir [3], the folder can be created to paste all species and cd is used to change the current folder [4]. With saveas, the graphs are saved in the required folder.
Program Output:
Plots for O2 species
1. Enthalpy vs Temperature
2. Entropy vs Temperature
3. Specific heat vs Temperature
Plots for N2 species
1. Enthalpy vs Temperature
2. Entropy vs Temperature
3. Specific heat vs Temperature
Plots for CO2 species
1. Enthalpy vs Temperature
2. Entropy vs Temperature
3. Specific heat vs Temperature
Screenshot of the window where all the folders are saved
Each folder out of 53 folders contains 3 plots (Enthalpy.jpg, Entropy.jpg and Specific heat.jpg).
A 'Molecular Weights.txt' file was created using fopen to save the molecular weight of each species calculated using the function 'molecularweight_calc.m'. The screenshot of 'Molecular Weights.txt' file is shown below:
Conclusion:
The theoretical concepts of thermodynamics and file parcing has been understood. A program in MATLAB has been created for file parcing. File Parsing is done successfully for the NASA Thermodynamic data and the thermodynamic properties such as enthalpy, entropy and specific heat has also found. The molecular weight of each species has also been found. The thermodynamic properties of each species has been successfully plotted with respect to temperature and is stored in a separate directory using MATLAB code.
References:
[1] https://www.computerhope.com/jargon/p/parse.htm
[2] https://in.mathworks.com/help/matlab/ref/str2num.html
[3] https://in.mathworks.com/help/matlab/ref/mkdir.html
[4] https://in.mathworks.com/help/matlab/ref/cd.html
[5] https://en.wikipedia.org/wiki/List_of_thermodynamic_properties#:~:text=Within%20thermodynamics%2C%20a%20physical%20property,which%20it%20can%20be%20specified.
[6] https://en.wikipedia.org/wiki/List_of_thermodynamic_properties#:~:text=Within%20thermodynamics%2C%20a%20physical%20property,which%20it%20can%20be%20specified.
[7] https://www.thermal-engineering.org/what-is-thermodynamic-property-definition/