Parsing NASA thermodynamic Data
ABSTRACT :
File parsing is very useful tool to manage a particular set of data use it for further studies and solving complex problems. In this project we have been provide NASA's thermodynamic data for various molecules (like O2, CO2, C2H6 etc) which contins values of coefficients of polynomial equations which are used to used to calculate specific heat, enthalpy and entropy. All these thermodynamic properties are function of temperature i.e. the value changes iven aas the temperature changes. We aim to parse the data given and plot the variation of each property as function of temperature for different molecules(also called as species).
Equations that were provided :
Cp `-=` Specific Heat
H `-=` Enthalpy
S `-=` Entropy
R `-=` Universal Gas Constant
ai `-=` Coefficients
T `-=` Temperature
Code for our problem :
1. Parsing the thermodynamic data to solve our problem
clear all
close all
clc
%Initial Parsing of data
f1 = fopen('THERMO.dat','r'); %opening file and get it in reading mode
fgetl(f1); %reading first line which is "THERMO'
temp_glob = fgetl(f1); %reading second line and storing the information in a variable
a = strsplit(temp_glob,' '); %spliting string into different parts and create a cell array of 1x4 size
%converting cell array elements into numbers
low_global_temp = str2double(a{2});
mid_global_temp = str2double(a{3});
high_global_temp = str2double(a{4});
%creating loop to read next 3 lines which start with "!" and not important for us
for i = 1:3
fgetl(f1);
end
%finding no. of species given in our data
%6th line contains 'O' which is 1st specie and after every 4 lines another
%specie name comes until line 214(sublime txt editor displays line numbers)
num_of_species = length((6:4:215));
%Creating for loop to plot graphs and calculate molecular weight for each species
for i = 1:num_of_species
specie_info = fgetl(f1);
b = strsplit(specie_info,' '); %finding the name of specie
specie_name = b{1};
p = strfind(specie_info,'.'); %finding local temperature
low_local_temp = str2double(specie_info(p(1)-3:p(1)+3));
high_local_temp = str2double(specie_info(p(2)-4:p(2)+3));
mid_local_temp = str2double(specie_info(p(3)-4:p(3)+3));
%reading all coefficent value form 1st line
coef_line_1 = fgetl(f1);
x = strfind(coef_line_1,'E');
a1 = str2double(coef_line_1(1:x(1)+3));
a2 = str2double(coef_line_1(x(1)+4:x(2)+3));
a3 = str2double(coef_line_1(x(2)+4:x(3)+3));
a4 = str2double(coef_line_1(x(3)+4:x(4)+3));
a5 = str2double(coef_line_1(x(4)+4:x(5)+3));
%reading all coefficent value form 2nd line
coef_line_2 = fgetl(f1);
y = strfind(coef_line_2,'E');
a6 = str2double(coef_line_2(1:y(1)+3));
a7 = str2double(coef_line_2(y(1)+4:y(2)+3));
a8 = str2double(coef_line_2(y(2)+4:y(3)+3));
a9 = str2double(coef_line_2(y(3)+4:y(4)+3));
a10 = str2double(coef_line_2(y(4)+4:y(5)+3));
%reading all coefficent value form 3rd line
coef_line_3 = fgetl(f1);
z = strfind(coef_line_3,'E');
a11 = str2double(coef_line_3(1:z(1)+3));
a12 = str2double(coef_line_3(z(1)+4:z(2)+3));
a13 = str2double(coef_line_3(z(2)+4:z(3)+3));
a14 = str2double(coef_line_3(z(3)+4:z(4)+3));
%Writing loop to calculate thermodyanmic properties for their temperature range
T = linspace(low_local_temp,high_local_temp,100);
R = 8.314;
for j = 1:length(T)
if T(j)<1000
cp = (a1 + a2.*T + a3.*T.^2 + a4.*T.^3 + a5.*T.^4).*R;
H = (a1 + (a2.*T)/2 + (a3.*T.^2)/3 + (a4.*T.^3)/4 + (a5.*T.^4)/5 + a6./T).*T.*R;
S = (a1*log(T)+ a2.*T+ (a3.*T.^2)/2 + (a4.*T.^3)/3 + (a5.*T.^4)/4 + a7).*R;
else cp = (a8 + a9.*T + 10.*T.^2 + a11.*T.^3 + a12.*T.^4).*R;
H = (a8 + a9.*T./2 + (a10.*T.^2)/3 + (a11.*T.^3)/4 + (a12.*T.^4)/5 + a13./T).*T.*R;
S = (a8*log(T)+ a9.*T+ (a10.*T.^2)/2 + (a11.*T.^3)/3 + (a12.*T.^4)/4 + a14).*R;
end
end
%plotting different graphs
%Ploting graph between Specific Heat(cp) and temperature
figure(1)
plot(T,cp,"Color",'b',"LineWidth",2);
xlabel('Temperature(K)')
ylabel('Specific Heat(cp)(J/Kg-K)')
title(specie_name)
%Ploting graph between Enthalpy and temperature
figure(2)
plot(T,H, "Color",'g','LineWidth',2);
xlabel('Temperature(K)')
ylabel('Enthalpy(H)(J/Kg)')
title(specie_name)
%Ploting graph between Entropy and temperature
figure(3)
plot(T,S,"Color",'m',"LineWidth",2);
xlabel('Temperature(K)')
ylabel('Entropy(J/K)')
title(specie_name)
%creating a folder and saving different plots in each one of them for
%different species
parentdir = '/MATLAB Drive/Species_Plots';
new_folder = fullfile(parentdir, sprintf('%s',specie_name));
status = mkdir(new_folder);
cd(new_folder)
saveas(1,'specific heat.jpg');
saveas(2,'Enthalpy.jpg');
saveas(3,'Entropy.jpg');
%generating molecular weight and specie name in command window
specie = specie_name
MM = molecular_weight(specie_name)
end
fclose(f1);Explaination of above code -
1. Under comment %Initial Parsing of data
Open our data file and read lines one by one using 'fgetl' function. Store global temperature range in suitable variables.
2. Next few lines of code are explained in comments itself.
3. Under comment %Creating for loop to plot graphs and calculate molecular weight for each species
First, split our current line using "strsplit" and save the name of molecule in variable named 'specie_name'. To save local temperature(unique for each specie), find the decimal point in the line using "strfind" and with suitable logic convert the string of temperature in numbers using "str2num". Next, we will read 3 lines containg coefficients and save all 14 coefficients in suitable variables.
4. Under comment %Writing loop to calculate thermodyanmic properties for their temperature range
We create a loop for evaluate value of each property at different temperature. Using if-else statement we determine value to use coefficients from a1 to a7 before 1000K temperature and use coeffecients from a8-a14 after 1000K. It must be noted that 1000K is common mid temperature for all species.
5. Under comment %plotting different graphs
Generate different graphs for each property.
6. Under comment %creating a folder and saving different plots in each one of them for different species
This part of program is written to save 3 plots for each specie in different folder with name of parent folder as 'Species_Plots'.
7. Under comment %generating molecular weight and specie name in command window
Calling the function named "molecular_weight" and passing name of specie as input .
2. Code for generating a function named "molecular_weight"
function [MM] = molecular_weight(specie_name)
st1 = specie_name;
MM = 0;
%Creating string arrays for basic elements involved in species
data_given = ['H','C','N','O','A','S'];
%Creating array of molecular weight of elements mentioned above at same respective
%position in the array
mole_wt = [1,12,14,16,40,32];
%creating for loop to calculate the molecular weight
for i = 1:length(st1)
for j = 1:length(data_given)
if strcmp(st1(i),data_given(j)) == 1
MM = MM + mole_wt(j);
u=j;
break;
end
end
if str2num(specie_name(i))>1
MM = MM + mole_wt(u)*(str2num(specie_name(i))-1);
end
endExplaination of above code -
1. Storing name of species in variable name 'st1'.
2. Creating string array with basic elements used in formation of species.
3. Creating array containg molecular weight of each element in previous array at same position.
4. Create loop to compare specie name with our string array generated inn point 2. and add one by one to calculate overall molecuular weight.
5. This step is necessary for species containg numbers and accounts that number in calculation of our molecular weight.
RESULTS:
1. Images(few) of Command Window displaying molecular weight of each species -
2. Image of all folders generated for each species -
3. Plots of O2 , N2 and CO2 -
i) Plots of O2
ii) Plots of N2
iii) Plots of CO2