Project - Parsing NASA thermodynamic data (MATLAB)
Aim: A MATLAB Code to parse the NASA thermodynamic data file and then calculate the thermodynamic properties of various gas species.
NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H, and S. They have also documented the coefficients that are required to evaluate these polynomials for 1000+ species.
- Write a function that extracts the 14 co-efficient and calculates the enthalpy, entropy, and specific heats for all the species in the data file. (Thermo.dat)
- Calculating the molecular weight of each species and showing in the command window.
- Plotting the Cp, Enthalpy, and Entropy for the local temperature range (low temperature: high temperature) specific for each species.
- Plots for O2, N2, and CO2
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, a5, a6, and a7 are the numerical coefficients provided in NASA thermodynamic files.
- The first 7 numbers starting on the second line of each species entry are the seven coefficients for the high-temperature range above 1000 K.
- The next seven numbers are the coefficients for the low-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 not equal to Delta Hf(T) and one needs to have the data for the reference elements to calculate Delta Hf(T).
- R is the universal gas constant = 314 J/(mol-K)
- T is the Temperature Range
Thus, for each species, there are 14 coefficients in all. In addition to the polynomial coefficients for each species, the Thermodynamic Data (Thermo.dat txt file) provides the species name, its elemental makeup, and the temperature ranges over which the fits are valid.
The important input to the fitting procedure is specific heat, enthalpy, and entropy as a function of temperature. In addition, the extent of the fit two temperature ranges, the temperature ranges have to be specified. Thus, the common temperature connecting the two ranges is 1000 K,
File Parsing: - File Parsing is a technique to read the information from a particular text file and that information will go to be used by MATLAB. If we know file parsing then we can read data from any other external file and add that data for our calculation. By doing we can perform data transformation. In other words, we can split large data into small data, and we can able to store them so that we can manipulate the data.
For this project, we are going to parsing NASA Thermodynamic format for CHEMKIN-II file which contains temperature ranges with 53 species. Thermo.dat file consists of sections containing element data, species data, thermo data, and reaction 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 the specified delimiter |
|
str2num |
Convert character array or string to numeric array |
|
strtok |
Return selected parts of the string |
|
strfind |
Find one string within another |
|
linspace |
Generate linearly spaced vectors |
|
fprintf |
Write formatted data to file |
|
pwd |
Identify current directory |
|
mkdir |
Make a new directory |
|
cd |
Change the working directory |
Steps Involve in Coding: -
- First of all, before coding the final NASA Thermodynamics, we have to create different functions related to Molecular weight, Read and writing molecular weight, Specific heat, Enthalpy, Entropy & Plotting.
- The function should be created by using if, else, and for loops.
- After creating the function code save the code with an appropriate name so that while calling the functions, we will not get any errors.
Detail Coding of functions
- Specific heat Calculation function code:
The specific heat is the amount of heat per unit mass required to raise the temperature by one degree 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;
end- Entropy 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;
end- Molecular Weight function code: -
Molecular weight is a measure of the sum of the atomic weight values of the atoms in a molecule. Molecular weight is used in chemistry to determine stoichiometry in chemical reactions and equations.
- Writing function for the molecular weight we have to specify a specific atomic number.
- The atomic number should match the atomic name.
- For this calculation here I am using
- Hydrogen = 1.008
- Carbon = 12.011
- Oxygen = 15.999
- Nitrogen = 14.007
- Argon = 39.948
- Now provide an initial starting range of molecular weight by 0.
- By using for loop add atomic number with each species to get particular 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 graphs 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 the program by clear and clean the workspace and command also close the background dialogue box.
- Next Step is to provide default directory by using ‘cd’ to start the program and save all the thing into that directory, here in the program I have used my default location.
- After that, we have to read the text file i.e Thermo.dat file and open it by using ‘fopen’ and r’’.
- Now for reading the line use ‘fgetl’.
Fig 1.1 Detail of Thermo.dat file
- Now next step is to obtain the temperature range from the Thermo.dat file.
- Now by using ‘strsplit’ we have to split the temperature into a particular string cell.
- After that, we have to use ‘str2num’ for converting string cell to numbers.
- Here we divided the temperature into a string and then converted to a number so that we can use separately to solve the NASA thermodynamics.
- Now in the text file, we have to skip the next three lines which we are not useful for coding. So, for that, we are using ‘for & end’ loops.
- After that, we have to read all the 53 species so that we can plot CP,H, and S. Also we can find the molecular weights.
- For reading all those 53 species here I am using ‘for’ loop.
- Now read the first line, here I use O_RL = fgetl.
- For reading specific species use ‘strtok’. This syntax will return to specific species after every reading.
- By ‘strfind’ we will use to find string find from another.
- Now we have to code for reading the local temperature range.
- By using linspace specify the temperature range from low to high and by using 1000 K.
- Also, add universal gas constant value ‘R’.
- Now next step is to extract higher and lower temperature coefficient by reading line by line from the species line.
- Now specify the specific heat, enthalpy, entropy for calling the function code.
- For calculating Molecular weight specify the code.
clear all; close all; clc
% Setting default folder on every startup of Matlab
cd 'D:\3D Design & Simulation\Skill lync\Challenges\Parsing NASA thermodynamic data\NASA Thermodynamic\Species'
...
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)
endOutput Results:
Plotting of O2 Graph
Fig 1.2 Specific Heat Vs Temperature Range
Fig 1.3 Enthalpy Vs Temperature Range
Fig 1.4 Entropy Vs Temperature Range
Plotting of N2 Graph
Fig 1.5 Specific Heat Vs Temperature Range
Fig 1.6 Enthalpy Vs Temperature Range
Fig 1.7 Entropy Vs Temperature Range
Plotting of Co2 Graph
Fig 1.8 Specific Heat Vs Temperature Range
Fig 1.9 Enthalpy Vs Temperature Range
Fig 1.10 Entropy Vs Temperature Range
- Screenshot of the folder location
- Molecular Weight Screenshot