Parsing NASA thermodynamic data
AIM - To write a MATLAB code to parse the NASA thermodynamic data and then to calculate
the thermodynamic properties for various gas species.
Objectives -
> To write a program that extracts 14 coefficients(first 7coefficients for high
temp & next 7 coefficients for low temp) and to calculate Specific heat, Enthalpy
& Entropy for each species in the data file.
> To calculate the molecular weight of all the species and display it in the
command window.
> To plot Cp, H & S for the local temperature range(low to high) specific for
each species.
> Save the plots as images with appropriate names and dump them in separate folders
for each species with a suitable name. All these should be generated
automatically.
File Parsing -
In general, 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. This aims at facilitating performing some kind of transformation of data. Parsing splits a file or other input into pieces of data that can be easily stored or manipulated.
NASA thermodynamic data -
NASA thermodynamic data i.e THERMO.dat file consists of 53 various gas species with the temperature range. Each species have 14 coefficients with a temperature range low to high. By parsing the data file of these 14 coefficients we calculate the thermodynamic properties i.e specific heat, enthalpy and entropy of each species.
Specific Heat(Cp) -
specific heat or heat capacity is a measure of the heat added to a mass of a system to produce a unit increase in temperature.
Cp = (a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4)*R
Enthalpy(H) -
Enthalpy is not a directly measurable property. It is a synthetic combination of the internal energy (U) and the flow energy (pV) exchanged with the surroundings.
H = (a1 + a2 T /2 + a3 T^2 /3 + a4 T^3 /4 + a5 T^4 /5 + a6/T)*R*T
Entropy(S) -
Entropy is a microscopic property associated with the microscopic energy transfer called heat (Q), between the system and its surroundings. Entropy is also called displacement of heat.It is not directly measurable, but can be related to other properties. Entropy is a measure of the level of irreversibility associated with any process.Unlike energy, it is nonconservative. It is a very important property in thermodynamics
S = (a1 lnT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7)*R
Where,
R = Universal Gas constant
T = Temperature
Cp = Specific heat[Kj]
H = Enthalpy[Kj/mol-K]
S = Entropy[Kj/mol]
a(i) = Coefficients[i = 1 to 14]
Program -
clear all
close all
clc
% Inputs
R = 8.314;
% Reading the first line of the file
f1 = fopen('Thermo.dat','r');
L1 = fgetl(f1); % Reading the first line of the file
% Reading the temperature of the file
L2 = fgetl(f1);
t = strsplit(L2); % Splitting
% Converting string cells to numbers
templow = str2num(t{2});
tempmid = str2num(t{3});
temphigh = str2num(t{4});
T = linspace(templow,temphigh,500); % Range of temperature
% Skipping next three lines
for i = 1:3
L1 = fgetl(f1);
end
% Reading the species
for i = 1:53
templine = fgetl(f1);
A = strsplit(templine,''); % Splitting into cell array
species = (A{1});
templine1 = fgetl(f1); % going to next line
a = strfind(templine1,'E'); % Finding the string 'E'
% Finding first 7 higher cofficients in first line
a1 = templine1(1:a(1)+3);
a1 = str2num(a1);
a2 = templine1(a(1)+4:a(2)+3);
a2 = str2num(a2);
a3 = templine1(a(2)+4:a(3)+3);
a3 = str2num(a3);
a4 = templine1(a(3)+4:a(4)+3);
a4 = str2num(a4);
a5 = templine1(a(4)+4:a(5)+3);
a5 = str2num(a5);
% Finding first 7 higher cofficients in second line
templine2 = fgetl(f1); % going to next line
b = strfind(templine1,'E'); % Finding the string 'E'
a6 = templine2(1:b(1)+3);
a6 = str2num(a6);
a7 = templine2(b(1)+4:b(2)+3);
a7 = str2num(a7);
% Finding second 7 lower cofficients in second line
a8 = templine2(b(2)+4:b(3)+3);
a8 = str2num(a8);
a9 = templine2(b(3)+4:b(4)+3);
a9 = str2num(a9);
a10 = templine2(b(4)+4:b(5)+3);
a10 = str2num(a10);
% Finding second 7 lower cofficients in third line
templine3 = fgetl(f1); % going to next line
c = strfind(templine1,'E'); % Finding the string 'E'
a11 = templine3(1:c(1)+3);
a11 = str2num(a11);
a12 = templine3(c(1)+4:c(2)+3);
a12 = str2num(a12);
a13 = templine3(c(2)+4:c(3)+3);
a13 = str2num(a13);
a14 = templine3(c(3)+4:c(4)+3);
a14 = str2num(a14);
mkdir(['/Users/nayanakumar/Documents/MATLAB/file parsing/',species])
% Calculate speafic heat, Enthalpy & Entropy
if (T<=tempmid)
Cp = R*(a8.*log(T)+(a9*T)+((a10*(T).^2)/2)+((a11*(T).^3)/3)+((a12*(T).^4)/4)+a14);
H = R*T.*(a8+((a9*T)/2)+((a10*(T).^2)/3)+((a11*(T).^3)/4)+((a12*(T).^4)/5)+((a13./T)));
S = R.*(a8.*log(T)+(a9.*T)+((a10*(T).^2)/2)+((a11*(T).^3)/3)+((a12*(T).^4)/4)+a14);
else
Cp = R*(a1.*log(T)+(a2*T)+((a3*(T).^2)/2)+((a4*(T).^3)/3)+((a5*(T).^4)/4)+a7);
H = R*T.*(a1+((a2*T)/2)+((a3*(T).^2)/3)+((a4*(T).^3)/4)+((a5*(T).^4)/5)+((a6./T)));
S = R.*(a1.*log(T)+(a2.*T)+((a3*(T).^2)/2)+((a4*(T).^3)/3)+((a5*(T).^4)/4)+a7);
end
cd(['/Users/nayanakumar/Documents/MATLAB/file parsing/',species])
%Plotting Cp v/s T
figure(1)
plot(T,Cp,'linewidth',3,'color','r');
xlabel('Temperature');
ylabel('Specific Heat');
title(sprintf('Specific Heat Vs Temperature range of %s',species));
grid on
% plotting Enthalpy Vs Temperature
figure(2)
plot(T,H,'linewidth',3,'color','b');
xlabel('Temperature');
ylabel('Enthalpy');
title(sprintf('Enthalpy Vs Temperature range of %s',species));
grid on
% Plotting Entropy Vs temperature
figure(3)
plot(T,S,'linewidth',3,'color','g');
xlabel('Temperature');
ylabel('Entropy');
title(sprintf('Entropy Vs Temperature range of %s',species));
grid on
saveas(1,'Cp_temp','jpg')
saveas(2,'H_temp','jpg')
saveas(3,'S_temp','jpg')
% Calculating the Molecular weight of each species
s_name = {'H' 'C' 'N' 'O' 'AR'}; % elements name
At_weight = [1 12 14 16 40]; % Atomic weight of elements
w = 0;
for i = 1:length(species)
for j = 1:length(s_name)
if strcmp(species(i),s_name(j))
w = w + At_weight(j);
position = j; % storing the position of matched element
end
end
n= str2num(species(i));
if n>1
w = w + At_weight(position)*(n-1)
end
end
fprintf('Molecular weight of %s :- ',species);
fprintf('%f',w);
disp(' ')
end
Explanation -
1. To begin with, we have given input of universal gas constant 'R'.
2. Then we have used 'fopen' command to open the thermo.dat file in the read mode i.e 'r' to read the data file.
3. To read the data line by line we use the command 'fgetl' , wherein the first line is the title/Head of the data file and we have ignored in the program.
4. Next line consists of Global temperature range, which is a string data hence we have split the data to get an individual range by using 'strsplit' command and convert split temperature range from string to number we have used 'str2num'.
5. Then we have given range for the temperature using 'linspace' command from temp low to temp high to the length of 500.
6. Then we have ignored the next three lines using 'for loop' command.
7. Then again 'for loop' is created for the 53 species in the data file and to calculate the said values of each species.
8. Next line in the 'for loop' we used 'fgetl' command to read the next line.
9. Once we read the line we again used the 'strsplit' command to split the data and then stored the first data in a file.
10. Then next three lines consists of coefficients of the high and low temperature.
11. To find each coefficient we used 'strfind' command and stored each coefficient in the numbers a1 to a14 also changed from string to number by using 'str2num' command.
12. That is high-temperature coefficients are from a1 to a7 and low-temperature coefficients are from a8 to a14.
13.In the next line, we created the folder to save all the species in each respective folder.
13. Then we have used 'if & else' command to calculate the Cp, H & S for the temperature range between low to high.
14. As we have calculated Cp, H & S for all species, the same have been plotted with respect to temperature and saved all the plots by creating the change in directory in their respective species folder.
15., At last, we calculated the Molecular weight of each species using 'for loop' command by giving input of element name and atomic weight of each element.
Reference -
1.http://combustion.berkeley.edu/gri-mech/data/nasa_plnm.html
2.http://www.thermix.net/2012/05/enthalpy-entropy-and-specific-heat.html
3.https://in.mathworks.com/matlabcentral/answers/140438-what-does-parse-error-basically-mean#:~:text=In%20general%2C%20the%20term%20%22parsing,when%20you%20execute%20the%20code.
4.https://en.wikipedia.org/wiki/Parsing
5.https://in.mathworks.com/matlabcentral/answers/106320-make-a-loop-to-calculate-molecular-weight