Project 1 - Parsing NASA thermodynamic data
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 & then calculate the thermodynamic properties of various gas speices.
Objective:
1. To write a function that extracts the 14 coefficients and calculates the enthalpy, entropy, and specific heats for all the species in the NASA thermodynamic data file "Thermo.dat".
The NASA polynomials have the form as shown below:
Where R is the universal gas constant, T is the temperature
2. To calculate the molecular weight of each species and display it in the command window.
3. To plot the Cp, Enthalpy, and Entropy for the local temperature range (low temperature: high temperature) specific for each species.
4. To save the plots as images with appropriate names and dump them in separate folders for each species with suitable name automatically.
File Parsing:
File parsing refers to how information can be read from a particular text file or any type of file and how that information can be used in a computer program.
Step 1 - Define the function
Step 2 - Create an Input Parser object
Step 3 - Add inputs to the scheme
Step 4 - Set properties to adjust parsing
Step 5 - Parse the inputs
Step 6 - Use the inputs in function
Step 7 - Calling function
NASA Thermodynamic File:
THERMO.dat file is a data file containing 14 co-efficients for 53 different elements given by NASA. These co-efficients are utilized to calculate the specific heat, enthalpy & entropy for a particular element. Each element has its own range of local temperature.
where,
R = Universal Gas constant
T = Temperature
Cp = Specific heat (KJ)
S = Entropy(KJ/mol-K)
H = Enthalpy(KJ/mol)
a(i) = Coefficients; where, i = 1 to 14
Specific Heat:
The specific heat is the amount of heat per unit mass required to raise the temperature by one degree celsius. The relationship between heat and temperature change is usually expressed in the form shown below where c is the specific heat.
Cp = R*(a1 + (a2*T) + (a3*(T).^2) + (a4*(T).^3) + (a5*(T).^4))
Entropy:
Entropy, the measure of a 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.
S = R*(a1.*log(T) + (a2*T) + ((a3*(T).^2)/2) + ((a4*(T).^3)/3) + ((a5*(T).^4)/4) + a7)
Enthalpy:
Enthalpy, a property of a thermodynamic system is equal to the system's internal energy plus the product of it's pressure and volume. In a system it's contained so as to prevent mass transfer, for processes at constant pressure, the heat absorbed or released equals the change in enthalpy.
H = R*T.*(a1 + ((a2*T)/2) + ((a3*(T).^2)/3) + ((a4*(T).^3)/4) + ((a5*(T).^4)/5) + ((a6./T)))
Main Code:
clear all
close all
clc
% Input
R = 8.314;
% Reading data file
f1 = fopen('THERMO.dat','r')
get = fgetl(f1);
temp = fgetl(f1);
s = strsplit(temp);
% Coverting strings into numbers
templow = str2num(s{2});
tempmid = str2num(s{3});
temphigh = str2num(s{4});
% Temperature range
T = linspace(templow,temphigh,1000);
% Skipping comments
for i=1:3
get = fgetl(f1);
end
% Reading the speices
for j = 1:53
tline = fgetl(f1);
% Splitting into cell array
A = strsplit(tline,' ');
% To read the species name
species = (A{1});
% Converting string to number
lowtemp = str2num(A{end-3});
midtemp = str2num(A{end-2});
hightemp = str2num(A{end-1});
tline1 = fgetl(f1)
a = strfind(tline1,'E');
% To find 1st 7 higher & 2nd 7 lower co-effs
% Finding higher co-effs
a1 = tline1(1:a(1)+3);
a1 = str2num(a1);
a2 = tline1(a(1)+4:a(2)+3);
a2 = str2num(a2);
a3 = tline1(a(2)+4:a(3)+3);
a3 = str2num(a3);
a4 = tline1(a(3)+4:a(4)+3);
a4 = str2num(a4);
a5 = tline1(a(4)+4:a(5)+3);
a5 = str2num(a5);
tline2 = fgetl(f1);
b = strfind(tline2,'E');
a6 = tline2(1:b(1)+3);
a6 = str2num(a6);
a7 = tline2(b(1)+4:b(2)+3);
a7 = str2num(a7);
% Finding lower co-effs
a8 = tline2(b(2)+4:b(3)+3);
a8 = str2num(a8);
a9 = tline2(b(3)+4:b(4)+3);
a9 = str2num(a9);
a10 = tline2(b(4)+4:b(5)+3);
a10 = str2num(a10);
tline3 = fgetl(f1);
c = strfind(tline3,'E');
a11 = tline3(1:c(1)+3);
a11 = str2num(a11);
a12 = tline3(c(1)+4:c(2)+3);
a12 = str2num(a12);
a13 = tline3(c(2)+4:c(3)+3);
a13 = str2num(a13);
a14 = tline3(c(3)+4:c(4)+3);
a14 = str2num(a14);
% Calculate specific heat
Cp = Specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,tempmid);
% Calculate entropy
S = Entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,tempmid);
% Calculate enthalpy
H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,R,tempmid);
% Calculate molecular weight of the speices
W = Molecular_weight(species);
% Writing molecular weight file
MW_file = fopen('Molecular Weight.txt','w');
fprintf(MW_file,'Molecular weight of %s is %d',species,W);
fclose(MW_file);
% Plotting Cp, H & S with Temperature
% Save it in a folder
cd('C:/Users/Student/Documents/MATLAB/species')
% Creating directory
mkdir(species);
% Changing directory
cd(species);
% Plot 1, Cp v/s T
figure(1)
plot(T,Cp,'linewidth',2,'color','b');
xlabel('Temperature(K)');
ylabel('Specific heat(KJ)');
grid on
title(sprintf('Specific heat v/s Temperature range of %s',species));
saveas(figure(1),'Specific Heat.jpg');
% Plot 2, S v/s T
figure(2)
plot(T,S,'linewidth',2,'color','r');
xlabel('Temperature');
ylabel('Entropy(KJ/mol-K)');
grid on
title(sprintf('Entropy v/s Temperature range of %s',species));
saveas(figure(2),'Entropy.jpg');
% Plot 3, H v/s T
figure(3)
plot(T,H,'linewidth',2,'color','g');
xlabel('Temperature');
ylabel('Enthalpy(KJ/mol)');
grid on
title(sprintf('Enthalpy v/s Temperature range of %s',species));
saveas(figure(3),'Enthalpy.jpg');
cd('C:/Users/Student/Documents/MATLAB')
% Molecular_weight
W = Molecular_weight(species)
endExplanation:
Firstly, we will Input the parameter R and we will open a file THERMO.dat file and reading this file using fopen command.
Then, we will read first THERMO using fgetl, strsplit command and then we read temperature file and we will split the temperatures using strsplit command to convert it into cell array and then we will convert this string to number using str2num command to find low, high and mid temperatures.
Then we will assign temperature range using linspace command.
Then we run the loop to skip command lines the data file using fgetl command inside the loop.
Now, we will run the loop to find 53 element in the data file by assigning fgetl command and then spliting into the cell array to find the species line and then find the species name and then find the number of the temperatures in the species line.
Then we read the new line. In order to read the co-efficient in the data file we will first find the 'E' string using strfind command to find before and after numbers of the co-efficients then we will find remaining 14 co-efficients and convert them to numbers from strings.
Out of which we will find first 7 Co-efficients which are for Higher temperature and then Second 7 Co-efficients which are for Lower temperature.
Now we will specify Specific heat Cp, Entropy S, Enthalpy H functions to find these parameters using these 14 co-efficients for Lower and Higher temperature
Then we will specify molecular weight W function to find the Molecular weight of the different species.
Now we will write the Molecular weight.txt file and fprintf of all the molecular weights of the species and then close the txt file using fclose command.
Now we change the directory using cd command and then creating directory folder where we have to save our files using mkdir command and then we will change the directory to this created folder.
Now, we will plot Specific heat, Entropy and Enthalpy with the Temperature for different species
We we will use scanf command in the title command to create different species title.
Then we will use saveas command to save the plots of different species in the directory folder.
Now we will change the directory to which our main code file runs in the MATLAB.
Function to find specific heat:
function Cp = Specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,tempmid)
if (T<= tempmid)
Cp = R*(a8 + (a9*T) + (a10*(T).^2) + (a11*(T).^3) + (a12*(T).^4));
else
Cp = R*(a1 + (a2*T) + (a3*(T).^2) + (a4*(T).^3) + (a5*(T).^4));
end
endNow we will find the Specific heat by creating Specific heat function and assigning it's related formula co-efficent from 14 co-efficients, gas constant, temperature and mid temperature.
Now, we will create a condition using if loop and make the temperature less than and equal to mid temperature to specify lower temperature co-efficent in the formula.
Then we will use else condition to specify co-efficient of higher temperature co-efficent in the formula.
Then end the loop and function to calculate Specific Heat.
Function to find Entropy:
function S = Entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,tempmid)
if (T <= tempmid)
S = R*(a8.*(log(T)) + (a9*T) + ((a10*((T).^2)./2)) + ((a11*((T).^3)./3)) + ((a12*((T).^4)./4)) +a14);
else
S = R*(a1.*(log(T)) + (a2*T) + ((a3*(T).^2)./2)) + ((a4*((T).^3)./3)) + ((a5*((T).^4)./4) + a7);
end
endNow we will find the Entropy by creating Entropy function and assigning its related formula coefficent from 14 coefficients, gas constant, temperature and mid temperature.
Now we will create a condition using if loop and make the temperature less than and equal to mid temperature to specify lower temperature co-efficent in the formula.
Then we will use else condition to specify coefficient of Higher temperature co-efficent in the formula.
Then end the loop and function to calculate Entropy.
Function to find Enthalpy :
function H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,R,tempmid)
if (T <= tempmid)
H = R*T.*(a8 + ((a9*T)/2)) + ((a10*((T).^2)./3)) + ((a11*((T).^3)./4)) + ((a12*(T).^4)./5)+ ((a13./T));
else
H = R*T.*(a1 + ((a2*T)/2) + ((a3*((T).^2)/3)) + ((a4*(T).^3)./4)) + ((a5*(T).^4)/5) + ((a6./T));
end
endNow we will find the Enthalpy by creating Enthalpy function and assigning it's related formula co-efficent from 14 co-efficients, Gas constant, Temperature and mid temperature.
Now we will create a condition using if loop and make the temperature less than and equal to mid temperature to specify lower temperature co-efficent in the formula.
Then we will use else condition to specify co-efficient of Higher temperature co-efficent in the formula.
Then end the loop and function to calculate Enthalpy.
Function to find Molecular weight of the different Species :
function W = Molecular_weight(species)
atomic_name = ['H','C','O','N','Ar'];
atomic_weight = [1.008 12.011 15.999 14.007 39.948];
W=0
for i=1:length(species)
for j=1:length(atomic_name)
if strcmp(species(i), atomic_name(j))
W = W + atomic_weight(j);
position = j;
end
end
n = str2num(species);
if n>1
W = W + atomic_weight(position)*(n-1);
end
end
fprintf('Molecular Weight of %s :- ',species);
fprintf('%f',W);
disp('')
endFirstly we create a Molecular_weight function W and assigning species input to that.
Then define the main element of all 53 species in an array atomic_name that are :- 'H','C','O','N','Ar'
Then by taking above main elements we will create there atomic_wieght array from any source like internet that are :-
1.008 12.011 15.999 14.007 39.948.
Then we will assign molecular starting range W To zero and then create a for loop condition which runs to the length of the species and an another nested for loop is created which runs to the length of the string atomic_name.
Then we will create a condition using if condition and use strcmp command to compare the species and atomic_name strings to calculate molecular masses if by comparison both string found same then we will assign its molecular weight.
Now, we will string to the number using str2num command assign to n.
Then we run a condition of if loop to find whether the species contains more than one element if there are more than 1 element that is n>1
Then we will increment the atomic weight by the multiplication of position and n-1 to the atomic_weight array.
Then we will print the Molecular weight using fprintf command and display it using disp command.
Screenshot of the window where all the folders are saved:
Output:
Plot for O2:
Plot for N2:
Plot for C02:
Output of Molecular weight of the Different Species:
Errors faced:
Main Code:
clear all
close all
clc
% Input
R = 8.314;
% Reading data file
f1 = fopen('THERMO.dat','r')
get = fgetl(f1);
temp = fgetl(f1);
s = strsplit(temp);
% Coverting strings into numbers
templow = str2num(s{2});
tempmid = str2num(s{3});
temphigh = str2num(s{4});
% Temperature range
T = linspace(templow,temphigh,1000);
% Skipping comments
for i=1:3
get = fgetl(f1);
end
% Reading the speices
for j = 1:53
tline = fgetl(f1);
% Splitting into cell array
A = strsplit(tline,' ');
% To read the species name
species = (A{1});
% Converting string to number
lowtemp = str2num(A{end-3});
midtemp = str2num(A{end-2});
hightemp = str2num(A{end-1});
tline1 = fgetl(f1)
a = strfind(tline1,'E');
% To find 1st 7 higher & 2nd 7 lower co-effs
% Finding higher co-effs
a1 = tline1(1:a(1)+3);
a1 = str2num(a1);
a2 = tline1(a(1)+4:a(2)+3);
a2 = str2num(a2);
a3 = tline1(a(2)+4:a(3)+3);
a3 = str2num(a3);
a4 = tline1(a(3)+4:a(4)+3);
a4 = str2num(a4);
a5 = tline1(a(4)+4:a(5)+3);
a5 = str2num(a5);
tline2 = fgetl(f1);
b = strfind(tline2,'E');
a6 = tline2(1:b(1)+3);
a6 = str2num(a6);
a7 = tline2(b(1)+4:b(2)+3);
a7 = str2num(a7);
% Finding lower co-effs
a8 = tline2(b(2)+4:b(3)+3);
a8 = str2num(a8);
a9 = tline2(b(3)+4:b(4)+3);
a9 = str2num(a9);
a10 = tline2(b(4)+4:b(5)+3);
a10 = str2num(a10);
tline3 = fgetl(f1);
c = strfind(tline3,'E');
a11 = tline3(1:c(1)+3);
a11 = str2num(a11);
a12 = tline3(c(1)+4:c(2)+3);
a12 = str2num(a12);
a13 = tline3(c(2)+4:c(3)+3);
a13 = str2num(a13);
a14 = tline3(c(3)+4:c(4)+3);
a14 = str2num(a14);
% Calculate specific heat
Cp = Specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,tempmid);
% Calculate entropy
S = Entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,tempmid);
% Calculate enthalpy
H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,R,tempmid);
% Calculate molecular weight of the speices
W = Molecular_weight(species);
% Writing molecular weight file
MW_file = fopen('Molecular Weight.txt','w');
fprintf(MW_file,'Molecular weight of %s is %d',species,W);
fclose(MW_file);
% Plotting Cp, H & S with Temperature
% Plot 1, Cp v/s T
figure(1)
plot(T,Cp,'linewidth',2,'color','b');
xlabel('Temperature(K)');
ylabel('Specific heat(KJ)');
grid on
legend('Specific heat v/s Temperature');
% Plot 2, S v/s T
figure(2)
plot(T,S,'linewidth',2,'color','r');
xlabel('Temperature');
ylabel('Entropy(KJ/mol-K');
grid on
legend('Entropy v/s Temperature');
% Plot 3, H v/s T
figure(3)
plot(T,H,'linewidth',2,'color','g');
xlabel('Temperature');
ylabel('Enthalpy(KJ/mol');
grid on
legend('Enthalpy v/s Temperature');
% Molecular_weight
W = Molecular_weight(species)
endFunction to find Specific heat
function Cp = Specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,tempmid)
if (T<= tempmid)
Cp = R*(a8 + (a9*T) + (a10*(T).^2) + (a11*(T).^3) + (a12*(T).^4));
else
Cp = R*(a1 + (a2*T) + (a3*(T).^2) + (a4*(T).^3) + (a5*(T).^4));
end
endFunction to find Entropy
function S = Entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,tempmid)
if (T <= tempmid)
S = R*(a8.*(log(T)) + (a9*T) + ((a10*((T).^2)./2)) + ((a11*((T).^3)./3)) + ((a12*((T).^4)./4)) +a14);
else
S = R*(a1.*(log(T)) + (a2*T) + ((a3*(T).^2)./2)) + ((a4*((T).^3)./3)) + ((a5*((T).^4)./4) + a7);
end
endFunction to find Enthalpy
function H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,R,tempmid)
if (T <= tempmid)
H = R*T.*(a8 + ((a9*T)/2)) + ((a10*((T).^2)./3)) + ((a11*((T).^3)./4)) + ((a12*(T).^4)./5)+ ((a13./T));
else
H = R*T.*(a1 + ((a2*T)/2) + ((a3*((T).^2)/3)) + ((a4*(T).^3)./4)) + ((a5*(T).^4)/5) + ((a6./T));
end
endFunction to find Molecular weight of different species
function W = Molecular_weight(species)
atomic_name = ['H','C','O','N','Ar'];
atomic_weight = [1.008 12.011 15.999 14.007 39.948];
W=0
for i=1:length(species)
for j=1:length(atomic_name)
if strcmp(species(i), atomic_name(j))
W = W + atomic_weight(j);
position = j;
end
end
n = str2num(species);
if n>1
W = W + atomic_weight(position)*(n-1);
end
end
fprintf('Molecular Weight of %s :- ',species);
fprintf('%f',W);
disp('')
endOutputs:
Due to the above code, the plots & the figures weren't correct
Couldn't figure out for which molecule the graphs had been plotted
The code wasn't saving the graphs in a folder in the same directory of the program files
So, I've made changes to the code & got the output correctly