Project 1 - Parsing NASA thermodynamic data
A MATLAB Program to Parsing NASA thermodynamic data
Objective: - A MATLAB Code to parse the NASA thermodynamic data file and then calculate thermodynamic properties of various gas species.
NASA came up with polynomials that can be used to test thermodynamic properties such as Cp, H and S. They also wrote the coefficients needed to analyze these correlations for the 1000+ variables.
Write a 14-function that interacts correctly and calculates enthalpy, entropy and specific temperatures for all the variables in the data file. (Thermo.dat)
The molecular weight of each species is calculated and displayed in the command window.
To plot the Cp, Enthalpy and Entropy of the local temperature range (low temperature: high temperature) specific to each type.
O2, N2 and CO2 sites
The formulae for calculating the NASA Thermodynamic are given:
`(Cp)/R = a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4`
`H/(RT) = a1 + (a2 T) /2 + (a3 T^2 )/3 + (a4 T^3 )/4 + (a5 T^4) /5 + (a6)/T`
`S/R = a1 lnT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7`
Where,
• a1, a2, a3, a4, a6, and a7 are numerical combinations provided by NASA's thermodynamic files.
• The first 7 digits starting in the second line of each entry are seven coefficients of the maximum temperature above 1000 K.
• The next seven digits are the coefficients of the temperature range below 1000 K.
• H in the above equation is defined as
H (T) = Delta Hf (298) + [H (T) - H (298)] so that, in general, H (T) is equal to Delta Hf (T) and one has to have data for Delta counters Hf (T).
• Gas is constant = 314 J / (mol-K)
• T is the temperature
Therefore, for each type there are 14 coefficients in all. In addition to the polynomial combinations of each type, Thermodynamic Data (Thermo.dat txt file) provides the name of the species, its basic composition, and the range of temperatures that are ready to work where applicable.
An important input to the proper process is the specific temperature, enthalpy, and entropy as a function of temperature. In addition, equilibrium causes two temperatures, temperatures must be specified. Therefore, the average temperature that connects the two distances is 1000 K,
File Parsing: - File Parsing is a way to read information from a specific file that such information will be used by MATLAB. If we know that a page entries in a file we can read the data from any other external file and add that data to our calculations. In practice we can make data changes. In other words, we can split big data into small data, and be able to store it so that we can use the information.
For this project we will transfer the NASA Thermodynamic format to the CHEMKIN-II file containing temperatures of 53 variants. The Thermo.dat file contains categories that contain object data, object types, thermo data, and response rate data.
Important Syntax Used during Coding
|
Fopen : |
Open a file or obtain information about open files |
|
‘r’ : |
Open the file for reading |
|
fgetl |
Read line from file |
|
strsplit |
Split string or character vector at specified delimiter |
|
str2num |
Convert character array or string to numeric array |
|
strtok |
Return selected parts of string |
|
strfind |
Find one string within another |
|
linspace |
Generate linearly spaced vectors |
|
fprintf |
Write formatted data to file |
|
pwd |
Identify current directory |
|
mkdir |
Make new directory |
|
cd |
Change working directory |
Steps Involve in Coding: -
• First, before writing the final NASA Thermodynamics, we have to create various functions related to Molecular Weight, Read and Write Molecular Weight, Direct Temperature, Enthalpy, Entropy and Plotting.
• Work should be performed using, once and for loops.
• After creating a work code save the code with the correct name so that while calling functions, we do not find any errors.
Detail Coding of functions
Specific heat Calculation function code:
The specific temperature is the amount of heat per unit required to raise the temperature by one Celsius.
`CP=(a1+a2*t+(a3*t^2)+(a4*t^3)+(a5*t^4))*R`
% Calculation of specific heat
function CP = Specificheat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp)
j = 1:length(Temp)
if Temp (j) <3000
end
% Comparing with global temperature
% If less than global it will take minimum coefficients
if Temp > global_medium_temp
CP =(a1+a2.*Temp+a3.*Temp.^2+a4.*Temp.^3+a5.*Temp.^4).*R;
% If temperature is more it will take maximum coefficients
else
CP =(a8+a9.*Temp+a10.*Temp.^2+a11.*Temp.^3+a12.*Temp.^4).*R;
end
Enthalpy function code:
A thermodynamic quantity is equivalent to the total heat content of a system. It is equal to the internal energy of the system plus the product of pressure and volume. Enthalpy is similar to energy, but not the same.
`H=(a8+(a9*t/2)+(a10*(t^2)/3)+(a11*(t^3)/4)+(a12*(t^4)/5)+ (a13/t))*R*t`
% Calculation of enthalpy
function H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,Temp,R,global_medium_temp)
j = 1:length(Temp)
if Temp (j) <3000
end
% Comparing with global temperature
% If less than global it will take minimum coefficients
if Temp > global_medium_temp
H =(a1.*Temp+(a2.*Temp.^2)./2+(a3.*Temp.^3)./3+(a4.*Temp.^4)./4+(a5.*Temp.^5)./5+a6).*R;
else
% If temperature is more it will take maximum coefficients
H =(a8.*Temp+(a9.*Temp.^2)./2+(a10.*Temp.^3)./3+(a11.*Temp.^4)./4+(a12.*Temp.^5)./5+a13).*R;
endEntropy function code:
A thermodynamic quantity representing the unavailability of a system's thermal energy for conversion into mechanical work often interpreted as the degree of disorder or randomness in the system. Entropy is also a measure of the number of possible arrangements the atoms in a system can have.
`s=(a8*log(t)+a9*t+(a10*(t^2)/2)+(a11*(t^3)/3)+(a12*(t^4)/4)+a14)*R`
% Calculation of entropy
function S = Entropy_calcEntropy_calc(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,Temp,R,global_medium_temp)
j = 1:length(Temp)
if Temp (j) <3000
end
% Comparing with global temperature
% If less than global it will take minimum coefficients
if Temp > global_medium_temp
S =(a1.*log(Temp)+a2.*Temp+a3.*(Temp.^2)./2+a4.*(Temp.^3)./3+a5.*(Temp.^4)./4+a7).*R;
else
% If temperature is more it will take maximum coefficients
S =(a8.*log(Temp)+a9.*Temp+a10.*(Temp.^2)./2+a11.*(Temp.^3)./3+a12.*(Temp.^4)./4+a14).*R;
endMolecular Weight function code:
Molecular weight is a measure of the total number of atomic weight in molecules. The molecular weight is used in chemistry to determine stoichiometry in chemical reactions and grade.
The molecular weighting function must specify a specific atomic number.
The atomic number must match the name of the atom.
For this figure here I use
Hydrogen = 1.008
Carbon = 12.011
Oxygen = 15.999
Nitrogen = 14.007
Argon = 39.948
Now provide the initial starting point for the molecular weight at 0.
Using the loop add the number of atoms and each species to obtain a specific molecular weight.
% Function to calculate molecular weight
function Molecularweight = Molecular_weight(Species_SP)
atomic_name = ['H','C','O','N','A']; % Standard Atomic Name
% Hydrogen, Carbon, Oxygen, Nitrogen, Argon
fprintf('Molecular Weight of ');
fprintf(Species_SP);
fprintf(' : ');
Atomic_Weight = [1.008 12.011 15.999 14.007 39.948]; % Standard Atomic Number
Molecularweight = 0; % Starting Range of Molecule
% Adding Atomic number to get particular value
for i = 1:length(Species_SP)
for j = 1:length(atomic_name)
if strcmpi(Species_SP(i),atomic_name(j))
Molecularweight = Molecularweight +Atomic_Weight(j);
position = j;
end
end
n = str2num(Species_SP(i));
if n>1
Molecularweight = Molecularweight + Atomic_Weight(position)*(n-1);
end
end
fprintf('molecular weight %s :',Species_SP);
fprintf('%d',Molecularweight);
disp(' ');
end
- Rear & Write molecular weight function
- This code is basically to create a text file for molecular weight.
% Read & Write Molecular text file
function Read_write (Species_SP,Molecularweight)
Molecular_list = fopen('Molecular weight.txt','a+');
fprintf( Molecular_list,'Molecular weight of %s is %d \n',Species_SP,Molecularweight);
fclose(Molecular_list);
end
- Plotting function code
- In plotting code, we will plot an individual CP, H, S Vs Temperature.
- Also, we will save all the graph to separate folders.
% Plotting for for each species
function plotting(CP,H,S,Temp,Species_SP)
% Saving as a folder
cur_dir=pwd;
mkdir(Species_SP);
cd(Species_SP);
% Plot 1, Specific heat (CP) vs Temperature Range (Temp)
figure(1)
plot(Temp,CP,'linewidth',2,'color','r');
xlabel('Temperature [K]');
ylabel('Specific Heat [kJ]');
title(sprintf('Specific Heat vs Temperature Range of %s',Species_SP));
grid on
saveas(figure(1),'Specific Heat.jpg');
% Plot 2, Enthalpy (H) vs Temperature Range (Temp)
figure(2)
plot(Temp,H,'linewidth',2,'color','b');
xlabel('Temperature [K]');
ylabel('Enthalpy in [kJ/(mol)]');
title(sprintf('Enthalpy vs Temperature Range of %s',Species_SP));
grid on
saveas(figure(2),'Ethalpy.jpg');
% Plot 3, Entropy (S) vs Temperature Range (Temp)
figure(3)
plot(Temp,S,'linewidth',2,'color','g');
xlabel('Temperature [K]');
ylabel('Entropy in [kJ/(mol-K)]');
title(sprintf('Entropy vs Temperature Range of %s',Species_SP));
grid on
saveas(figure(3),'Entropy.jpg');
cd(cur_dir);
% Molecular weight
M = Molecular_weight(Species_SP)
f2 = fopen('molecular mass','w');
fprintf('%s n %s n %f n','species','Molecular mass',M)
fclose(f2)
end
- Steps Involved in Main Code:
Start writing a program by deleting and cleaning the workspace and ordering and closing the domain dialog box.
The Next Step is to provide the default directories by using 'cd' to start the program and save everything to that directory, here in the program I used my default location.
After that, we should read the text file i.e Thermo.dat file and open it using 'fopen' and r ''.
Now for reading the line use 'fgetl'.
Now the next step is to find the temperature range in the Thermo.dat file.
Now by using 'strsplit' we have to divide the temperature into a cell of a particular cell.
After that, we'll have to use 'str2num' to convert a character cell to numbers.
Here we split the temperatures into rope and converted to a number so that we can use them separately to solve NASA thermodynamics.
Now in the text file, you must skip the next three rows that are not useful for encoding. So, with that, we use the 'end & end' traps.
After that, we have to study all 53 variables so that we can program CP, H, and S. And we can get the molecular weights.
By reading all those 53 genres here I am using the 'for' loop.
Now read the first line, here I use O_RL = fgetl.
To learn some types of words 'strtok'. This syntax will return to certain types after reading it all.
With 'strfind' we will use to find the thread for another.
Now we will have to write about reading the local temperature.
Using linspace describe the temperature range from low to high and use 1000 K.
Also, add the amount of permanent gas 'R'.
Now the next step is to extract the maximum and minimum temperature by reading line by line from the species line.
Now describe the specific heat, the enthalpy, the entropy of the cost function code.
Molecular weight calculator define the code.
clear all; close all; clc
% Setting default folder on every startup of Matlab
cd 'C:\Users\Rockstar Sudip\Desktop\e'
...
pwd
% Reading Data file - NASA Thermodynamic
f1 = fopen ('THERMO.dat','r'); % Open a file & 'r' Open the file for reading
gt = fgetl(f1);
% Obtaing Temperatures (Reading Header Information)
global_temp = fgetl (f1);
temp_range = ('300.000 1000.000 5000.000'); % Spliting
A = strsplit(temp_range,' ');
% Now convert string cells tonumbers using str2num,
global_low_temp = str2num (A{1}); % Global Low
global_medium_temp = str2num (A{2}); % Global Medium
global_high_temp = str2num (A{3}); % Global High
% skipping the next three comments in the file
for i=3:5
comments = fgetl(f1);
end
for i = 1:53
O_RL = fgetl(f1); % RL = Read line
Species_SP = strtok(O_RL); % SP = Split
L_T = strfind (O_RL,'G');
% reading local temperatue range
Avg_temp = O_RL(L_T+1:end-1);
Avg_temp = str2num(Avg_temp);
% Range of temperature & Gas Constant
Temp = linspace(global_low_temp,global_high_temp,1000); % Temperature Range assumed
R = 8.314; % Universal Gas Constant J/(mol-K)
% Extracting higher and lower temperature coeffs
line_1 = fgetl(f1);
line_2 = fgetl(f1);
line_3 = fgetl(f1);
A = strfind(line_1,'E');
B = strfind(line_2,'E');
C = strfind(line_3,'E');
% Extracting HIGH-temperature coefficients - FIRST 7 coefficients
a1 = str2num(line_1(1:A(1)+3));
a2 = str2num(line_1(A(1)+4:A(2)+3));
a3 = str2num(line_1(A(2)+4:A(3)+3));
a4 = str2num(line_1(A(3)+4:A(4)+3));
a5 = str2num(line_1(A(4)+4:A(5)+3));
a6 = str2num(line_2(1:B(1)+3));
a7 = str2num(line_2(B(1)+4:B(2)+3));
% Extracting LOW-temperature coefficients - SECOND 7 coefficients
a8 = str2num(line_2(B(2)+4:B(3)+3));
a9 = str2num(line_2(B(3)+4:B(4)+3));
a10 = str2num(line_2(B(4)+4:B(5)+3));
a11 = str2num(line_3(1:C(1)+3));
a12 = str2num(line_3(C(1)+4:C(2)+3));
a13 = str2num(line_3(C(2)+4:C(3)+3));
a14 = str2num(line_3(C(3)+4:C(4)+3));
% Specific heat calculation
CP = Specificheat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp);
% Enthalpy calculation
H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,Temp,R,global_medium_temp);
% Entropy calculation
S = Entropy_calc(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,Temp,R,global_medium_temp);
% Calculating molecular weight of species
Molecularweight(i) = Molecular_weight(Species_SP);
% Read & Write Molecular list
Read_write(Species_SP,Molecularweight)
% Creating individual folder for each species and plotting
Plotting(CP,H,S,Temp,Species_SP)
end
Output Results:
Plotting of O2 Graph
Plotting of N2 Graph
Plotting of Co2 Graph
Screenshot of the folder location:
Molecular Weight Screenshot
Reference
- http://www.thphys.nuim.ie/CompPhysics/matlab/help/techdoc/ref/fopen.html
- http://matlab.izmiran.ru/help/techdoc/ref/fgetl.html
- http://matlab.izmiran.ru/help/techdoc/ref/strtok.html
- http://matlab.izmiran.ru/help/techdoc/ref/strtok.html
- https://www.mathworks.com/help/matlab/ref/strsplit.html
- https://www.mathworks.com/help/matlab/ref/str2num.html
- http://matlab.izmiran.ru/help/techdoc/ref/strfind.html
- http://matlab.izmiran.ru/help/techdoc/ref/pwd.html