NASA thermodynamic file parsing
Aim : To write a function that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the data file (NASA thermodynamic data file provided in the project)
Project : We will be using concept of fileparsing in matlab to open and read the file, then extract the coefficients and finally plot data sets for different species.
Explanation of the code:
clear all
close all
clc
%opening file and extracting header information
f1 =fopen(\'THERMO.dat\',\'r\')
first_line = fgetl(f1)
second_line = fgetl(f1)
%extracting data from second line
A = strsplit(second_line,\' \')
global_min_temp = str2num(A{2})
global_mid_temp = str2num(A{3})
global_high_temp = str2num(A{4})
%ignoring comments
for i=1:3
fgetl(f1)
end
In this section of the code ,all we did was to open the file named THERMO.dat usng fopen command. We used a file handle f1. After that we used fgetl command to read the lines.Now the second line contained global temperature values which we needed to extract. To do that, we used strsplit (string split commmand) and double spaces as splitting logic to extract the vaues as string. The string values needed to be converted into number by using str2num (string to number) command. Next three lines were comments, which were extracted by using simple for loop.Then , we proceeded to the species block in next section of the code.
line1 = fgetl(f1)
species = strsplit(line1,\' \')
species_name(k) = species(1)
temp = line1(49:74)
B = strsplit(temp,\' \')
local_low_temp = str2num(B{1})
local_high_temp = str2num(B{2})
local_mid_temp = str2num(B{3})
T = linspace(local_low_temp,local_high_temp,200);
tline = fgetl(f1)
a = findstr(tline,\'E\')
a1= tline(1:a(1)+3);
a1 = str2num(a1)
a2 = tline(a(1)+4:a(2)+3)
a2 = str2num(a2)
a3 = tline(a(2)+4:a(3)+3)
a3 = str2num(a3)
a4 = tline(a(3)+4:a(4)+3)
a4 = str2num(a4)
a5 = tline(a(4)+4:a(5)+3)
a5 = str2num(a5)
tline2 = fgetl(f1)
a6 = tline2(1:a(1)+3)
a6 = str2num(a6)
a7 = tline2(a(1)+4:a(2)+3)
a7 = str2num(a7)
a8 = tline2(a(2)+4:a(3)+3)
a8 = str2num(a8)
a9 = tline2(a(3)+4:a(4)+3)
a9 = str2num(a9)
a10 = tline2(a(4)+4:a(5)+3)
a10 = str2num(a10)
tline3 = fgetl(f1)
a11 = tline2(1:a(1)+3)
a11= str2num(a11)
a12 = tline2(a(1)+4:a(2)+3)
a12 = str2num(a12)
a13 = tline2(a(2)+4:a(3)+3)
a13 = str2num(a13)
a14 = tline2(a(3)+4:a(4)+3)
a14 = str2num(a14)
This part of the code deals with the extraction of high temperature and low temperature coefficients from the species block. There were , in total, 53 different species and the entire process was repeated 53 times using for loop. First 7 coefficients were high temperature coefficients and next 7 were low temperature coefficients.
Now we parsed next line using fgetl command again. This line consisted of species name and local temperature ranges, so we used strsplit command and double spaces as logic to extract the name of the species. The temperatures were present at positions ranging from 49:74 in the line, so we first extracted these and again used strsplit command and then str2num command to convert string values of temperature to numerical values. The temperatures were saved as local_low_temp , local_mid_temp and local_high_temp. After that, we divided the temperature values into 200 values ranging from low to high using linspace command.
After this, we moved on to the coefficients. Next line was read again using fgetl command. Now we used the logic that next coefficeint was 3 places away from \'E\' string in the code. So we used findstr command. \'E\' string was at 12,27,32,47,62nd place in each of the lines; each species consisting of 3 lines, line 1 and 2 containing 5 constants each and line3 containing 4 elements. So using this logic, we extracted all the 14 coefficents of the species. After extraction of the constants , we moved to the calculation part.
These are the formulas we used to calculate speciic heat(Cp), Enthalpy(H), Entropy(S) and R is the universal gas constant
for j =1:length(T)
%high temp values
if T>local_mid_temp
Cp(j) = (a1 + a2*T(j) + a3*T(j).^2 + a4*T(j).^3 +a5*T(j).^4)*R;
H(j) = (a1 + (a2*T(j))/2 + (a3*T(j)^2)/3 + (a4*T(j)^3)/4 +(a5*T(j)^4)/5 + (a6/T(j)))*R*T(j)
S(j) = (a1*log(T(j))+a2*T(j)+(a3*T(j)^2)/2+(a4*T(j)^3)/3+ (a5*T(j)^4)/4+a7)*R
%low temp values
else
Cp(j) = (a8 + a9*T(j) + a10*T(j).^2 + a11*T(j).^3 +a12*T(j).^4)*R;
H(j) = (a8 + (a9*T(j))/2 + (a10*T(j).^2)/3 + (a11*T(j).^3)/4 +(a12*T(j).^4)/5 + (a13/T(j)))*R*T(j)
S(j) = (a8*log(T(j))+a9*T(j)+(a10*T(j).^2)/2+(a11*T(j).^3)/3+ (a12*T(j).^4)/4+a14)*R
end
end
Using a for loop and if else statements, we calculated the values. We calculated values for low temperature and high temperature separately; usng local_mid_temp as reference value.After calulcation we moved onto creating separate folders and then plotting the results and storing them in the creating folders.
dataname = species_name(k)
sp_name = char(dataname)
mkdir([\'C:\\Users\\abhishek\\Documents\\nasa\\\',sp_name]);
cd([\'C:\\Users\\abhishek\\Documents\\nasa\\\',sp_name]);
figure(1)
title [dataname]
plot(T,Cp)
xlabel(\'Temperature(in K)\')
ylabel(\'Specific_heat(in J/Kg-K)\')
saveas(gca,\'temperature vs specific heat\', \'jpeg\');
cd(pwd)
figure(2)
title [dataname]
plot(T,H)
xlabel(\'Temperature(in K)\')
ylabel(\'enthalpy(J/Kg)\')
saveas(gca,\'temperature vs enthalpy\', \'jpeg\');
cd(pwd)
figure(3)
title [dataname]
plot(T,S)
xlabel(\'Temperature(in K)\')
ylabel(\'entropy(in J/Kg-K)\')
saveas(gca,\'temperature vs entropy\', \'jpeg\');
cd(pwd)
We had calculated the species name earlier, species(k) for each loop, this was saved in another variable dataname. Then we used char command to get the characters from dataname and saved them as sp_name which will be the names of our folders. Then we used the commands mkdir and cd to create and change the directories in which our results will be stored.
FINAL CODE
clear all
close all
clc
%opening file and extracting header information
f1 =fopen(\'THERMO.dat\',\'r\')
first_line = fgetl(f1)
second_line = fgetl(f1)
%extracting data from second line
A = strsplit(second_line,\' \')
global_min_temp = str2num(A{2})
global_mid_temp = str2num(A{3})
global_high_temp = str2num(A{4})
%ignoring comments
for i=1:3
fgetl(f1)
end
%parsing species block
for k = 1:53
line1 = fgetl(f1)
species = strsplit(line1,\' \')
species_name(k) = species(1)
temp = line1(49:74)
B = strsplit(temp,\' \')
local_low_temp = str2num(B{1})
local_high_temp = str2num(B{2})
local_mid_temp = str2num(B{3})
T = linspace(local_low_temp,local_high_temp,200);
tline = fgetl(f1)
a = findstr(tline,\'E\')
a1= tline(1:a(1)+3);
a1 = str2num(a1)
a2 = tline(a(1)+4:a(2)+3)
a2 = str2num(a2)
a3 = tline(a(2)+4:a(3)+3)
a3 = str2num(a3)
a4 = tline(a(3)+4:a(4)+3)
a4 = str2num(a4)
a5 = tline(a(4)+4:a(5)+3)
a5 = str2num(a5)
tline2 = fgetl(f1)
a6 = tline2(1:a(1)+3)
a6 = str2num(a6)
a7 = tline2(a(1)+4:a(2)+3)
a7 = str2num(a7)
a8 = tline2(a(2)+4:a(3)+3)
a8 = str2num(a8)
a9 = tline2(a(3)+4:a(4)+3)
a9 = str2num(a9)
a10 = tline2(a(4)+4:a(5)+3)
a10 = str2num(a10)
tline3 = fgetl(f1)
a11 = tline2(1:a(1)+3)
a11= str2num(a11)
a12 = tline2(a(1)+4:a(2)+3)
a12 = str2num(a12)
a13 = tline2(a(2)+4:a(3)+3)
a13 = str2num(a13)
a14 = tline2(a(3)+4:a(4)+3)
a14 = str2num(a14)
R = 8.314
for j =1:length(T)
%high temp values
if T>local_mid_temp
Cp(j) = (a1 + a2*T(j) + a3*T(j).^2 + a4*T(j).^3 +a5*T(j).^4)*R;
H(j) = (a1 + (a2*T(j))/2 + (a3*T(j)^2)/3 + (a4*T(j)^3)/4 +(a5*T(j)^4)/5 + (a6/T(j)))*R*T(j)
S(j) = (a1*log(T(j))+a2*T(j)+(a3*T(j)^2)/2+(a4*T(j)^3)/3+ (a5*T(j)^4)/4+a7)*R
%low temp values
else
Cp(j) = (a8 + a9*T(j) + a10*T(j).^2 + a11*T(j).^3 +a12*T(j).^4)*R;
H(j) = (a8 + (a9*T(j))/2 + (a10*T(j).^2)/3 + (a11*T(j).^3)/4 +(a12*T(j).^4)/5 + (a13/T(j)))*R*T(j)
S(j) = (a8*log(T(j))+a9*T(j)+(a10*T(j).^2)/2+(a11*T(j).^3)/3+ (a12*T(j).^4)/4+a14)*R
end
end
dataname = species_name(k)
sp_name = char(dataname)
mkdir([\'C:\\Users\\abhishek\\Documents\\nasa\\\',sp_name]);
cd([\'C:\\Users\\abhishek\\Documents\\nasa\\\',sp_name]);
figure(1)
title [dataname]
plot(T,Cp)
xlabel(\'Temperature(in K)\')
ylabel(\'Specific_heat(in J/Kg-K)\')
saveas(gca,\'temperature vs specific heat\', \'jpeg\');
cd(pwd)
figure(2)
title [dataname]
plot(T,H)
xlabel(\'Temperature(in K)\')
ylabel(\'enthalpy(J/Kg)\')
saveas(gca,\'temperature vs enthalpy\', \'jpeg\');
cd(pwd)
figure(3)
title [dataname]
plot(T,S)
xlabel(\'Temperature(in K)\')
ylabel(\'entropy(in J/Kg-K)\')
saveas(gca,\'temperature vs entropy\', \'jpeg\');
cd(pwd)
end
Results: 1) Plots for O2
2) Plots for N2
3) Plots for CO2
Dump folder
Errors faced:
1) Problems were faced while saving folders because of incorrect snytax
Materials referred:
1) mathsworks.com