Parsing NASA thermodynamic data in MATLAB
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. The formulae are shown below.
Image 1.1 Formulae to calculate the enthalpy, entropy and specific heat.
The first 7 coefficients are HIGH-temperature coefficients and the SECOND 7 coefficients are LOW-temperature coefficients.
Also to use the local temperature range which is unique for each individual species, not the global temperature range.
1. To calculate the molecular weight of each species and display it in the command window.
2. To plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species.
3. To save the plots as images with appropriate names and dump them in separate folders for each species with suitable name.
4. Also the main program should take the "THERMO.dat" file as an input and generate Cp, H and S plots.
5. To write a detailed report to describe the logic used.
THEORY: About the given THERMO.dat file
- The given data file is thermochemical data file, used with another data file named grimech30.dat, used in the Optimization mechanism designed to model natural gas combustion, developed by NASA, named GRI-Mech 3.0.
- GRI-Mech 3.0 is an optimized mechanism designed to model natural gas combustion, including NO formation and reburn chemistry.
- The optimization process is designed to provide sound basic kinetics which also furnish the best combined modeling predictability of basic combustion properties.
- Improvements were made in the categories of updating the kinetics with recent literature results, including some new and improved target experiments to the optimization, expanding the mechanism and target selection, and examining the sensitivity to the thermodynamics.
The data file format is as shown below
PROCEDURE:
Process Flow:
- The program that is to be developed should do following operations
- Read the data from the THERMO.dat file.
- Store the Global high, low and mid temperature.
- Store the Species name as character and calculate its molecular mass.
- Store its local high, mid and low temperature.
- Store the 14 thermal coefficients.
- Calculate the Specific heat, enthalpy and Entropy for each and every elements or Compounds.
- Plot the Temperature vs Specific heat, enthalpy and Entropy graph and store all the graphs in separate folder for every species.
- Separate variables are used to store the global temperatures.
The flow of the program is shown below
- Clear all the stored data before.
- Open the file using the 'fopen' syntax.
- Storing the Global low, mid and high temperature by splitting the line and using the location number to store the temperature.
- Reading the species name by splitting the string and saving the name for further use.
- In that same line the local temperatures are present and is stored in some variable.
- Storing the thermal coefficients by storing the value of the coefficient using the location of the number.
- Creating 3 function for the calculation of specific heat, enthalpy and entropy.
- Calling the function to calculate the above.
- Creating a new folder which store the image of graph
- Plotting the graph of the temperature vs specific heat, enthalpy and entropy.
- Storing the graph as image.
- Calculating the molecular mass.
- Printing the value of the molecular mass in command window.
Step 1: Clearing the all data stored before in workspace and command window.
Below syntax is used to clear the data stored in workspace and the command window.
clear all
close all
clc- Clear all is to clear all the data stored in the workspace.
- Close all is to close all the figures opened.
- clc is used to clear the command window.
Step 2: Opening the file THERMO.dat file
Below syntax is used to open the THERMO.dat file and store the file handle in a variable.
f1 = fopen('THERMO.dat','r');In the above syntax 'f1' stores the file handle and opens the THERMO.dat file and 'r' is used for Read mode i.e it only read the file.
Step 3: Storing the Global temperatures.
Below syntax is used to store the Global low, mid and high temperature.
fgetl(f1);
FL1 = fgetl(f1);
FLS = strsplit(FL1);
global_low_temp = str2num(FLS{2});
global_mid_temp = str2num(FLS{3});
global_high_temp = str2num(FLS{4});- 'fgetl(f1)' is a syntax used to go to next line.
- Now FL1 store the next line and it is split using "strsplit".
- The three variables store in the syntax stores the three temperatures.
- Below image shows the image of the data file where the global temperatures are mentioned.
Image 1.2 THERMO.dat file
Step 4: Reading the species name by splitting the string and saving the name for further use.
Below syntax shows the syntax used to read the species name by splitting the line and saving the name in a variable.
line_001 = fgetl(f1);
a = strsplit(line_001);
Species_Name = (a{1});In the above syntax line_001 store the 6th line of the data file and a store the string from the line using the "strsplit".
The first string is the name of the element or compound and hence it is stored as Species_name.
Step 5: Storing the local low, mid and high temperature.
Below syntax shows the syntax used to store the local temperature for each compound or element.
ltmin = str2num(a{end-3});
ltmax = str2num(a{end-2});
ltmid = str2num(a{end-1});
The 3 variables store the local low, mid and high temperature for every compound or element.
Below image shows the image of the data of local temperature and species name.
Image 1.3 Image of the line in data file
Step 6: Storing the thermal coefficients by storing the value of the coefficient using the location of the number.
Below syntax is used to store the thermal coefficients in the form of an array.
%Extraction of the thermal coefficients in the first line
x = 1;
line_002 = fgetl(f1);
a = findstr(line_002, 'E');
for j = 1:length(a);
b = line_002(a(j)-11:a(j)+3);
F(x) = str2num(b);
x = x + 1;
end
%Extraction of the thermal coefficients in the second line
line_003 = fgetl(f1);
a = findstr(line_003, 'E');
for j = 1:length(a)
b = line_003(a(j)-11:a(j)+3);
F(x) = str2num(b);
x = x+1;
end
%Extraction of the thermal coefficients in the third line
line_004 = fgetl(f1);
a = findstr(line_004, 'E');
for j = 1:length(a)
b = line_004(a(j)-11:a(j)+3);
F(x) = str2num(b);
x = x+1;
end- There were 14 thermal coefficients displayed in 3 lines, and in first and second line there were 5 coefficients and in the last line there was only 4 coefficients.
- The "for" loop is used to extract the thermal coefficients in each line. Therfeore there were 3 loops in total to extract the coefficients in three lines.
- First the line is stored in line_002 and the position of E is found using the "findstr" syntax.
- Then the using the location of the "E" the 3 characters from the after "E" and 11 characters before "E" is stored in an Array named "F".
- Initially a variable named "x" is defined and is used as the location of the coefficient in the array and added 1 to it when it stores the coefficient so that next coefficient is stored in x+1 location of F.
- Below image shows the first line of the coefficients in the data file.
Image 1.4 Thermal coefficients in the data file.
Step 7: Creating 3 function for the calculation of specific heat, enthalpy and entropy.
Instead of putting the equation to calculate the specific heat, enthalpy and entropy in the same program function is created to simplify the program.
Calculation of Specific Heat(Cp)
Below equation shows the formula to calculate the Specific heat
`Cunderset(p)( ) = R*(aunderset(1)( ) + aunderset(2)( )*T + aunderset(3)( ) * T^2 + aunderset(4)( )* T^3 + aunderset(5)( )*T^4)`
Below syntax shows the function to calculate the Specific Heat using the thermal coefficients.
function [CP] = cp(a1,a2,a3,a4,a5,T)
R = 8.314;
CP = R*(a1 + (a2*T) + a3*(T^2) + a4*(T^3) + (a5*(T^4)));
end- In the above function named "cp" receives input a1,a2,a3,a4,a5 and T where a1 to a5 are the thermal coefficients and T is the temperature and gives the specific heat as output.
- Universal gas constant is also defined in the function.
Calculation of Enthalpy(H)
Below equation shows the formula to calculate the Enthalpy
`H = R*T*(aunderset(1)( ) + aunderset(2)( )*T /2+ aunderset(3)( ) * T^2/3 + aunderset(4)( )* T^3/4 + aunderset(5)( )*T^4/5 + a underset(6)( )*T^-1)`
Below syntax shows the function to calculate the Enthalpy(H) using the thermal coefficients.
function [H] = enthalpy(a1,a2,a3,a4,a5,a6,T)
R = 8.314;
H = R*T*(a1 + (a2*T/2) + (a3*(T^2)/3) + (a4*(T^3)/4) + (a5*(T^4)/5) + a6/T);
end- In the above function named "enthalpy" receives input a1,a2,a3,a4,a5,a6 and T where a1 to a6 are the thermal coefficients and T is the temperature and gives the enthalpy as output.
Calculation of Entropy(S)
Below equation shows the formula to calculate the Entropy
`S = R*(aunderset(1)( ) lnT + aunderset(2)( )*T + aunderset(3)( ) * T^2/2 + aunderset(4)( )* T^3/3 + aunderset(5)( )*T^4/4 + a underset(7)( ))`
Below syntax shows the function to calculate the Entropy(S) using the thermal coefficients.
function [S] = entropy(a1,a2,a3,a4,a5,a7,T)
R = 8.314;
lnT = log(T);
S = R*(a1*lnT + a2*T + (a3* T^2)/2 + (a4*T^3)/3 + (a5*T^4)/4 + a7);
end- In the above function named "entropy" receives input a1,a2,a3,a4,a5,a7 and T where a1 to a7 are the thermal coefficients and T is the temperature and gives the entropy as output.
Step 8: Calling the function to calculate the above
Below syntax is used to calculate the specific heat, enthalpy and entropy for the temperature range i.e form local low temperature to high temperature.
a1 = F(1);
a2 = F(2);
a3 = F(3);
a4 = F(4);
a5 = F(5);
a6 = F(6);
a7 = F(7);
a8 = F(8);
a9 = F(9);
a10 = F(10);
a11 = F(11);
a12 = F(12);
a13 = F(13);
a14 = F(14);
temperature_range = [ltmin:ltmax];
for j = 1:length(temperature_range)
Entropy(j) = 0;
Enthalpy(j) = 0;
Specific_heat(j) = 0;
if temperature_range(j) >= ltmid
Specific_heat(j) = cp(a1,a2,a3,a4,a5,temperature_range(j));
Enthalpy(j) = enthalpy(a1,a2,a3,a4,a5,a6,temperature_range(j));
Entropy(j) = entropy(a1,a2,a3,a4,a5,a7,temperature_range(j));
else
Specific_heat(j) = cp(a8,a9,a10,a11,a12,temperature_range(j));
Enthalpy(j) = enthalpy(a8,a9,a10,a11,a12,a13,temperature_range(j));
Entropy(j) = entropy(a8,a9,a10,a11,a12,a14,temperature_range(j));
end
end- Initially the all the thermal coefficient is temporarily stored in the a1,a2...a14 respectively.
- The first 7 coefficients are called as high temperature coefficients and next 7 are low temperature coefficients.
- A for loop is used to calculate the enthalpy, entrophy and specific heat for entire temperature range.
- 'if-else' condition is used as the high temperature coefficients are to be used if the temperature is above the local mid temperature.
- For the entire temperature range the results of it is stored in the form of array.
Step 9: Creating a new folder which store the image of graph
Below syntax is used to create the folder which is further used to store the images of the graph.
mkdir(Species_Name);- The "mkdir" is the command used to create the folder in the current location of the program.
- This syntax creates a folder in the name of the element or compound which is stored in the Species_Name.
Step 10: Plotting the graph of the temperature vs specific heat, enthalpy and entropy.
Below syntax is used to plot the graph of Temperature vs Specific heat save it as image in the folder created.
% Plotting the Temperature vs Specific heat
plot(temperature_range,Specific_heat)
url = sprintf('E:/Subramanya G S Skill-Lync/MATLAB_for_ME/Project_001/%s/Specific Heat Plot.jpg',Species_Name);
title(Species_Name)
xlabel('Temperature[K]')
ylabel('Specific Heat[kJ/kg]')
saveas(gcf,url)
% Plotting the Temperature vs Enthalpy
plot(temperature_range,Enthalpy)
url1 = sprintf('E:/Subramanya G S Skill-Lync/MATLAB_for_ME/Project_001/%s/Enthalpy Plot.jpg',Species_Name);
title(Species_Name)
xlabel('Temperature[K]')
ylabel('Enthalpy[kJ/mol]')
saveas(gcf,url1)
%plotting the Temperature vs Entropy
plot(temperature_range,Entropy)
url2 = sprintf('E:/Subramanya G S Skill-Lync/MATLAB_for_ME/Project_001/%s/Entropy Plot.jpg',Species_Name);
title(Species_Name)
xlabel('Temperature[K]')
ylabel('Entropy[J/K]')
saveas(gcf,url2)- Plot command plots the graph.
- sprintf is the command used format data into string or character vector
- The title, xlabel and ylabel is added in the program to display the x-axis and yaxis and the title for the graph.
Step 11: Storing the graph as image.
Below syntax is used to save the graph in image format.
saveas(gcf,url)- "saveas" saves the current frame using the "gcf" command and stores in the folder, whose location is saved in "url", in ".jpj" format.
Step 12: Calculating the molecular mass.
Below syntax is used to calculate the molecular mass of the element or compound stored in the Species_Name.
tline = Species_Name;
X = char('C','H','O','N','A');
Z = [12,1,16,14,40];
N = char('1','2','3','4','5','6','7','8','9');
c = tline == X ;
n = tline == N;
n1 = find(n == 1);
c1 = find(c == 1);
M = 0;
for j = 1:width(c)
for l = 1:length(X)
if find(c(l,j)==1) == 1
M = M + Z(l);
end
end
end
for l = 1:length(n1)
n2 = mod(n1(l),9);
n3 = (n1(l)-n2)/9;
c2 = find(c(:,n3) == 1) ;
M = M + (n2-1)*Z(c2);
end- The variable tline stores the Species_Name which had stored the element or compound name mentioned in the data file.
- The elements used in the data file is carbon(C), oxygen(O), Hydrogen(H), Nitrogen(N) and Argon(Ar).
- The characters of the elements used in the compound in the data file is stored in the form of characters array.
- In another array the molecular mass of the element is stored.
- And in the last array the numbers ranging from 1 to 9 is stored in the form of characters.
- "c" is a logical array that stores 1 if the condititon satisfies and 0 if not, i.e stores 1 if the the characters in tline belong to array X and 0 if not.
- Similarly "n" stores 1 if the character in the "tline" belong the to the array N.
- And n1 stores the location of 1 in the logical array "n", c1 stores the location of the 1 in the logical array "c".
- Initially M, which stores the molecular mass, is cleared.
- In the first for loop the mass is calculated only for the characters i.e if there is compound "CH4" the mass is calculated is only for CH and the 4th power the H is not calculated.
- This power is calculated in the second for loop.
- The first loop find the 1's in the logical array "c" and adds the number in the 'x' location of "Z", where x is equal to the row of logical array "c" in which the 1 is found.
- In the second loop "n1" gives the value of the power for the element i.e for a compound "CH4", n1 = 4.
- And "n3" gives the location of the element whose mass is to be added 'n' times, where n is power of element in the compound.
- Using the "n3" and comparing with the logical array "c" the element is found.
- Now the mass of the element which has n power in compound is add (n-1) time to the mass calculated in first loop.
- The mass is added n-1 times because its mass is added once in the first loop.
Step 13: Printing the value of the molecular mass in command window.
Below syntax is used to print the molecular mass of the compound
format = 'The molecular mass of %s is %d';
sprintf(format,Species_Name ,M)- In the above syntax the "The molecular mass of %s is %d" is the sentence which is stored in "format" where %s and %d is space that can replaced by character and number in respective place.
- sprintf prints the "format" and the %s is replaced with current Species_Name and %d is replaced with M value.
Step 14: Running the above procedure for 53 time
- This program should run for 53 time as there are 53 elements or compound in teh data file.
- Hence for loop is used to run this program for 53 time.
PROGRAM: The complete porgram is shown below
%Program to read the data in the given NASA GRI-Mech Version 3.0
%Thermodynamics file.
%clearing the data stored in workspace
close all
clear all
clc
%Opening the file THERmOD
f1 = fopen('THERMO.dat','r'); % Opening the THERMO.dat file
fgetl(f1); % Skipping first line
FL1 = fgetl(f1); % Storing 2nd line of the data file
FLS = strsplit(FL1); % splitting the same line
global_low_temp = str2num(FLS{2}); % Converting the string to number and storing the same
global_mid_temp = str2num(FLS{3});
global_high_temp = str2num(FLS{4});
% Initial conditions
F = 0;
x = 1;
% Skipping the 3 lines of comments in the data file
fgetl(f1);
fgetl(f1);
fgetl(f1);
% Starting the for loop to store the data from the file and compute the
% functions according to the challenge description
for i = 1:53
line_001 = fgetl(f1); % Storing the First line of the actual data of the species
a = strsplit(line_001); % Splitting the line
Species_Name = (a{1}); % The first part of the string is the name of the Species and is saved
ltmin = str2num(a{end-3}); % The 4th part of the string is the local low temperature
ltmax = str2num(a{end-2}); % The 5th part of the string is the local high temperature
ltmid = str2num(a{end-1}); % The 6th part of the string is the local mid temperature
%Storing the first five coefficients
line_002 = fgetl(f1); % Storing the second line of the actual data of the species
a = findstr(line_002, 'E'); % Finding the position of the "E" in the data
% For loop to store the first five coefficients stored in the second line of the actual data
for j = 1:length(a);
b = line_002(a(j)-11:a(j)+3); % Storing the character form 12 characters before "E" till three characters after "E"
F(x) = str2num(b); % Converting the string to number and storing the same in the "x"th location of array "F"
x = x + 1; % Adding 1 to it so that the next coefficient will be stored in the "x+1" location in array "F"
end
% Storing the second five coefficients
line_003 = fgetl(f1);
a = findstr(line_003, 'E');
% For loop to store the first five coefficients stored in the second line of the actual data
for j = 1:length(a)
b = line_003(a(j)-11:a(j)+3); % Storing the character form 12 characters before "E" till three characters after "E"
F(x) = str2num(b); % Converting the string to number and storing the same in the "x"th location of array "F"
x = x+1; % Adding 1 to it so that the next coefficient will be stored in the "x+1" location in array "F"
end
%Storing the last five coefficients
line_004 = fgetl(f1);
a = findstr(line_004, 'E');
for j = 1:length(a)
b = line_004(a(j)-11:a(j)+3); % Storing the character form 12 characters before "E" till three characters after "E"
F(x) = str2num(b); % Converting the string to number and storing the same in the "x"th location of array "F"
x = x+1; % Adding 1 to it so that the next coefficient will be stored in the "x+1" location in array "F"
end
x = 1; % "x" value reaches 14 therefore it needs to be defined again for the next set of data for different species
% Storing the coefficients in separate variable for further calculation
% of the Specific heat, Enthalpy and entropy
a1 = F(1);
a2 = F(2);
a3 = F(3);
a4 = F(4);
a5 = F(5);
a6 = F(6);
a7 = F(7);
a8 = F(8);
a9 = F(9);
a10 = F(10);
a11 = F(11);
a12 = F(12);
a13 = F(13);
a14 = F(14);
% Storing the temperature ranging from low to high temperature to low
% temperature
temperature_range = [ltmin:ltmax];
for j = 1:length(temperature_range)
Entropy(j) = 0;
Enthalpy(j) = 0;
Specific_heat(j) = 0;
if temperature_range(j) >= ltmid % if contition to use the high temperature coefficients if the temperature is more than the mid temperature
Specific_heat(j) = cp(a1,a2,a3,a4,a5,temperature_range(j)); % Calling the "cp" function to calculate the specific heat
Enthalpy(j) = enthalpy(a1,a2,a3,a4,a5,a6,temperature_range(j)); % Calling the "enthalpy" function to calculate enthalpy
Entropy(j) = entropy(a1,a2,a3,a4,a5,a7,temperature_range(j)); % Calling the "entropy" function to calculate entropy
else
Specific_heat(j) = cp(a8,a9,a10,a11,a12,temperature_range(j));
Enthalpy(j) = enthalpy(a8,a9,a10,a11,a12,a13,temperature_range(j));
Entropy(j) = entropy(a8,a9,a10,a11,a12,a14,temperature_range(j));
end
end
% Plotting the graph and storing the same as jpj image format
mkdir(Species_Name); % Creating the folder in the name of the characters stored in the Species_Name
% Plotting the Temperature vs Specific Heat
plot(temperature_range,Specific_heat)
% Storing the file location and file name
url = sprintf('E:/Subramanya G S Skill-Lync/MATLAB_for_ME/Project_001/%s/Specific Heat Plot.jpg',Species_Name);
title(Species_Name) % Storing the file location and file name
xlabel('Temperature[K]') % Label for X axis
ylabel('Specific Heat[kJ/kg]') % Label for y axis
saveas(gcf,url) % Saving the current frame and storing the same in the previously defined folder location
% Plotting the Temperature vs Enthalpy
plot(temperature_range,Enthalpy)
% Storing the file location and file name
url1 = sprintf('E:/Subramanya G S Skill-Lync/MATLAB_for_ME/Project_001/%s/Enthalpy Plot.jpg',Species_Name);
title(Species_Name) % Storing the file location and file name
xlabel('Temperature[K]') % Label for X axis
ylabel('Enthalpy[kJ/mol]') % Label for y axis
saveas(gcf,url1) % Saving the current frame and storing the same in the previously defined folder location
% Plotting the Temperature vs Entropy
plot(temperature_range,Entropy)
% Storing the file location and file name
url2 = sprintf('E:/Subramanya G S Skill-Lync/MATLAB_for_ME/Project_001/%s/Entropy Plot.jpg',Species_Name);
title(Species_Name) % Storing the file location and file name
xlabel('Temperature[K]') % Label for X axis
ylabel('Entropy[J/K]') % Label for y axis
saveas(gcf,url2) % Saving the current frame and storing the same in the previously defined folder location
% Clearing the data stored from the above calculations for next set of
% data
Specific_heat = 0;
Entropy = 0;
Enthalpy = 0;
% Calculation of the Molecular Mass of the element or compound
tline = Species_Name; % Storing the character of Species_Name in the tline
X = char('C','H','O','N','A'); % Storing the Characters of the elements symbol
Z = [12,1,16,14,40]; % Storing the elements mass in respective position as that of the "X"
N = char('1','2','3','4','5','6','7','8','9'); % Storing the numbers in the form of characters
c = tline == X; % Comparing the tline with X for finding the elements symbol in the chemical
n = tline == N; % Comparing the tline with N for finding the any power for the element
n1 = find(n == 1); % Finding the location of 1 in the logical array "c"
c1 = find(c == 1); % Finding the location of 1 in teh logocal array "n"
% Initial condition, Molecular mass = M
M = 0;
% For loop for calculating the Molecular mass irrespective of their powers
for j = 1:width(c)
for l = 1:length(X)
if find(c(l,j)==1) == 1 % if condititon which find 1 in "c"
M = M + Z(l); % Adds the number in the "l"th location of array Z
end
end
end
% for loop to add the mass of the element which has powers to it
for l = 1:length(n1)
n2 = mod(n1(l),9); % Finds the power
n3 = (n1(l)-n2)/9; % Finds the location of element which had power
c2 = find(c(:,n3) == 1) ; % Finds location of that element in the array Z
M = M + (n2-1)*Z(c2); % Calculates the total mass
end
% printing the Molecular mass of the element or compound
format = 'The molecular mass of %s is %d'; % Format of statement to print in the command window
sprintf(format,Species_Name ,M) % Prints the Statement
end
RESULTS:
On running this program
- it read the THERMO.dat file
- Stored the data in the respective variables
- Calculated specific heat, enthalpy and entropy
- the molecular mass of the each species is displayed in the command window
- and there were 53 sets of 3 graphs, namely Temperature vs Specific heat, Temperature vs Enthalpy and Temperature vs Entropy, stored in 53 different folder with name as chemical formula for which graph are plotted
CONCLUSION: A program to do the below instructions is written and explaination to the code is also mentioned above.
- to read the data from the given THERMO.dat file,
- to calculate the enthalpy, specific heat and entropy,
- to plot the graphs and store them in image format and
- to display the molecular mass in the command window.
The program has successfully gave the required results and all the graphs are attached below.