Parsing NASA thermodynamic data using Matlab Program
Aim: To write a code in Matlab to parse the NASA Thermodynamic data file and then calculate the Thermodynamic properties of various gas species.
Objective:
1. Write the function that extracts the 14 co-efficients and calculates the Enthalpy,Entropy and Specific heats of all species in the NASA thermodynamic data "THERMO.dat" file.
The NASA Polynomials have the form as shown below:
Where R is the universal gas constant, T is the temperature
The First 7 coefficients are High-temperature co-efficients and the Second 7 coefficients are Low-temperature coefficients.
2. Calculate the molecular weight of each species and display it in the command window.
3. Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature:high temperature) specific for each species.
4. Save the plots as images with appropriate names and dump them in seperate folders for each species with suitable name.All these should be generated automatically.
Background: NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S. They have also documented the co-efficients that are required to evaluate these polynomials for 1000+ species.
The thermodynamic properties in terms of the co-efficients are given as:
Where,
Cp - Specific heat at constant pressure (J/Kg-k)
H - Enthalpy (Joules - J)
S - Entropy (J/K)
T - Temperature (Kelvin)
R - Universal gas constant (R=8.314 J/mol-k)
(a1,a2,.....,a7) - Coefficients of species
Definitions:
Thermodynamics: Thermodynamics is a branch of physics that deals with heat, work, and temperature, and their relation to energy, radiation, and physical properties of matter.The behavior of these quantities is governed by the four laws of thermodynamics which convey a quantitative description using measurable macroscopic physical quantities, but may be explained in terms of microscopic constituents by statistical mechanics. Thermodynamics applies to a wide variety of topics in science and engineering, especially physical chemistry, chemical engineering and mechanical engineering, but also in other complex fields such as metrology.
Specific Heat:The specific heat capacity,(Cp) of a substance is the heat capacity of a sample of the substance divided by the mass of the sample. Informally, it is the amount of energy that must be added, in the form of heat, to one unit of mass of the substance in order to cause an increase of one unit in its temperature.The SI unit of specific heat is joule per kelvin and kilogram, J/(K kg).
Enthalpy: Enthalpy is a property of a thermodynamic system, that is a convenient state function preferred in many measurements in chemical, biological, and physical systems at a constant pressure. It is defined as the sum of the system's internal energy and the product of its pressure and volume.
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. The concept of entropy provides deep insight into the direction of spontaneous change for many everyday phenomena.
Molecular Mass: The molar mass is defined as the mass of a given substance divided by the amount of a substance and is expressed in g/mol.
Molecular Weight: The definition of molecular weight is most authoritatively synonymous with molecular mass; however, in common practice, it is also highly variable as are the units used in conjunction with it. Many common preparatory sources use g/mol and effectively define it as a synonym of molar mass, while more authoritative sources use Da or u and align its definition more closely with the molecular mass.
File Parsing: File parsing refers to how information cab be read from a particular text file or any type of file and how that information can be used in a computer program.
NASA Thermodynamic file: THERMO.dat file is a data file containing 14 coefficients for 53 different elements given by NASA. These coefficients are utilized to calculate specific heat, enthalpy & entropy for a particular element. Each element has its own range of local temperature.
List of 53 Elements contained in 'THERMO.dat' file:
AR C C2H C2H2 C2H3 C2H4 C2H5 C2H6 C3H7 C3H8 CH CH2
CH2(S) CH2CHO CH2CO CH20 CH2OH CH3 CH3CHO CH3O CH3OH CH4 CN CO
CO2 H H2 H2CN H2O H2O2 HCCO CCOH HCN HCNN HCNO HCO
HNCO HNO HO2 HOCN N N2 N20 NCO NH NH2 NH3 NNH
NO NO2 O O2 OH
CODE FOR PARSING THE DATA AND EVALUATE MOLECULAR WEIGHT FOR SPECIES:
MAIN CODE:
%%*** PARSING THE THERMODYNAMIC DATA AND EVALUATE MOLECULAR WEIGHT FOR SPECIES ***
% Clearing all previous memory in workspace
clear all
% Terminating all the previous plots
close all
% Clearing all the values in command window
clc
% Importing the 'THERMO.dat' file contents
ntd = fopen('THERMO.dat','r');
Header = fgetl(ntd);
% Reading Global Temperature Line
Global_temp = fgetl(ntd);
c = strsplit(Global_temp,' ');
global_low_temp = str2double(c{2});
global_mid_temp = str2double(c{3});
global_high_temp = str2double(c{4});
% Skipping 3 lines
fgetl(ntd);
fgetl(ntd);
fgetl(ntd);
% Loop for all 53 elements
for i = 1:53
line = fgetl(ntd);
% Starting with the first species
B = strsplit(line,' ');
% Obtaining the species name
species_name = (B{1});
% Finding the location of 'G' string
L = strfind(line,'G');
% Obtaining the local temperature values
local_low_temp = line(L+3:55);
local_high_temp = line(L+13:65);
local_mid_temp = line(L+23:75);
lt = str2double(local_low_temp);
ht = str2double(local_high_temp);
mt = str2double(local_mid_temp);
% Extracting the 14 co-efficients
% Line1 - Reading the coefficient data
coeffline1 = fgetl(ntd);
% Finding the location of 'E' string
X = strfind(coeffline1,'E');
% Extracting all the co-efficients from 1st line
a1 = coeffline1(1:X(1)+3);
a1 = str2double(a1);
a2 = coeffline1(X(1)+4:X(2)+3);
a2 = str2double(a2);
a3 = coeffline1(X(2)+4:X(3)+3);
a3 = str2double(a3);
a4 = coeffline1(X(3)+4:X(4)+3);
a4 = str2double(a4);
a5 = coeffline1(X(4)+4:X(5)+3);
a5 = str2double(a5);
% Line2 - Reading the coefficient data
coeffline2 = fgetl(ntd);
% Finding the location of 'E' string
Y = strfind(coeffline2,'E');
% Extracting all the co-efficients from 2nd line
a6 = coeffline2(1:Y(1)+3);
a6 = str2double(a6);
a7 = coeffline2(Y(1)+4:Y(2)+3);
a7 = str2double(a7);
a8 = coeffline2(Y(2)+4:Y(3)+3);
a8 = str2double(a8);
a9 = coeffline2(Y(3)+4:Y(4)+3);
a9 = str2double(a9);
a10 = coeffline2(Y(4)+4:Y(5)+3);
a10= str2double(a10);
% Line3 - Reading the coefficient data
coeffline3 = fgetl(ntd);
% Finding the location of 'E' string
Z = strfind(coeffline3,'E');
% Extracting all the co-efficients from 3rd line
a11 = coeffline3(1:Z(1)+3);
a11 = str2double(a11);
a12 = coeffline3(Z(1)+4:Z(2)+3);
a12 = str2double(a12);
a13 = coeffline3(Z(2)+4:Z(3)+3);
a13 = str2double(a13);
a14 = coeffline3(Z(3)+4:Z(4)+3);
a14 = str2double(a14);
% Local temperature range of the species
T = linspace(lt,ht,250);
% Specific Heat
Cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,mt,T);
% Entropy
S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,mt,T);
% Enthalpy
H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,mt,T);
NASA=sprintf('%s',species_name);
% Assigning the current directory to a variables 'cur_dir'
cur_dir = pwd;
% To creat folders for species
mkdir(NASA);
cd(NASA);
% Plotting Graphs
% Plot 1. Specific Heat vs Temperature
figure(1);
plot(T,Cp,'linewidth',2,'color','r')
xlabel('Temperature in K');
ylabel('Specific Heat in KJ');
title('Specific Heat vs Temperature');
grid on
% Plot 2. Enthalpy vs Temperature
figure(2);
plot(T,H,'linewidth',2,'color','b')
xlabel('Temperature in K');
ylabel('Enthalpy in KJ/(mol)');
title('Enthalpy vs Temperature');
grid on
% Plot 3. Entropy vs Temperature
figure(3);
plot(T,S,'linewidth',2,'color','g')
xlabel('Temperature in K');
ylabel('Entropy in KJ/(mol K)');
title('Entropy vs Temperature');
grid on
% Saving the plots
saveas(figure(1),'Cp.png');
saveas(figure(2),'H.png');
saveas(figure(3),'S.png');
cd(cur_dir);
% For Calculating Molecular Weight of the species
species_elements = {'H','C','N','O','A'};
AtomicWeight = [1 12 14 16 40];
M = 0;
for k = 1:length(species_name)
for j = 1:length(species_elements)
if strcmp(species_name(k),species_elements(j))==1
M = M+AtomicWeight(j);
p = 1;
end
end
b = str2double(species_name(k));
% Adding the weight of the repeating species_elements
if b>1
M = M+(AtomicWeight(p))*(b-1);
end
end
% Printing the calculated molecular weights
fprintf('Molecular Weight of [%s] - %damu n',species_name,M)
end
Program Process:
1. Compiler start the program
2. Then the data file is being imported to Matlab ("THERMO.dat" file).
3. Reading global temperature line.
4. Next the first five lines are being skipped.
5. The loop start for all the 53 species.
6. The 6th line of data file is being read.
7. The line is being split by ' ' (Space) delimiter.
8. The species name is stored in variable called 'species_name'
9. Finding the location of 'G' string in the line.
10. Obtaining the local temperature values
11. Extracting all the 14 coefficients and reading the 1st coefficient data line.
12. Finding the 'E' string in the line and extracting all the co-efficients from 1st line.
13. Reading the coefficients of line 2 and finding 'E' string then extracting all the co-efficients from 2nd line.
14. Reading the coefficients of line 3 and finding 'E' string then extracting all the co-efficients from 3rd line.
15. Functions for Specific heat, enthapy and entropy.
16. Assigning the variables to the current directory 'cur_dir' and creating the folder for species.
17. Then the graphs are plotted Cp,H & S vs Specific Temperature.
18. The plots are saved in the individual folders.
19. Next the molecular weight is calculated for each species.
20. The compiler terminates the program.
Syntax Used during coding:
Fopen: Open a file or obtain information about open files.
'r': Open the file for reading.
fgetl: Reading line from line.
strsplit: Split string or character vector at specified delimiter.
str2num: Convert character array or string to numeric array.
str2double: Convert string to double precision value.
linspace: Generate linearly spaced vectors.
strfind: Find one string within another
fprintf: Write formatted data to file.
pwd: Identify current directory.
mkdir: Make new directory.
cd: Change new directory.
CODE FOR CALCULATING THERMODYNAMIC PROPERTIES:
FUNCTION CODE: SPECIFIC HEAT(Cp)
function Cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,mt,T)
R = 8.314;
for i = 1:length(T)
if T(i)>=mt
Cp(i) = (a1+(a2*T(i))+(a3*(T(i)^2))+(a4*(T(i)^3))+(a5*(T(i)^4)))*R;
else
Cp(i) = (a8+(a9*T(i))+(a10*(T(i)^2))+(a11*(T(i)^3))+(a12*(T(i)^4)))*R;
end
endFUNCTION CODE: ENTROPY(S)
function S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,mt,T)
R = 8.314;
for i = 1:length(T)
if T(i)>=mt
S(i) = ((a1*log(T(i)))+(a2*T(i))+((a3*(T(i)^2))/2)+((a4*(T(i)^3))/3)+((a5*(T(i)^4))/4)+a7)*R;
else
S(i) = ((a8*log(T(i)))+(a9*T(i))+((a10*(T(i)^2))/2)+((a11*(T(i)^3))/3)+((a12*(T(i)^4))/4)+a14)*R;
end
endFUNCTION CODE: ENTHALPY(H)
function H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,mt,T)
R = 8.314;
for i = 1:length(T)
if T(i)>=mt
H(i) = R*T(i)*((a1)+(a2*(T(i))/2)+((a3*(T(i)^2))/3)+((a4*(T(i)^3))/4)+((a5*(T(i)^4))/5)+(a6/T(i)));
else
H(i) = R*T(i)*((a8)+(a9*(T(i))/2)+((a10*(T(i)^2))/3)+((a11*(T(i)^3))/4)+((a12*(T(i)^4))/5)+(a13/T(i)));
end
endProgram Process for all Functions:
1. The compiler start the program.
2. It checks the all input parameters from the main code.
3. Then according to the conditional statements the values of Cp,H & S are calculated.
4. The complier terminates the program.
OUTPUT:
The plots for all 53 species were generated successfully and saved in their respective folders.
O2 Plots:
N2 Plots:
Co2 Plots:
Plots Folder Window:
Molecular Weight:
CONCLUSION:
The NASA data file was imported and the thermodynamic properties were calculated successfully. Then plot the graphs for Cp,H&S vs Specific temperature and then calculated the molecular weight for each species.