Challenge

Project 1 - Parsing NASA thermodynamic data

AAkshat Hardaha15 Jul 2022

Parsing NASA thermodynamic data
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. 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 InputParser 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 coefficients for 53 different elements given by
NASA. These coefficients are utilized to calculate the specific heat, enthalpy & entropy for a
particular element. Each element has its own range of local temperature. Below screenshot
of the THERMO.dat file is displayed.

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
CH30
CH3OH
CH4
CN
CO
C02
H
H2
H2CN
H20
H202
HCCO
CCOH
HCNN
HCNO
HCO
HNCO
HNO
H02
HOCN
N
N2
N20
NCO
NH
NH2
NH3
NNH
NO
NO2
02
OH

Enthalpy, Entropy, and Specific Heat: Enthalpy, entropy, and specific heat, of various species values at different temperatures are required extensively to analyze the thermodynamic process.
For example, if the combustion process is considered during this, these quantities are required to observe how the combustion is going on and how efficient it will be. It eventually allows seeking the scope of improvement in the process. Similarly in the other processes like refrigeration, air conditioning, chemical kinetics, nozzle flow, detonation process, etc, these quantities are widely involved to observe the phenomenon. NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H, and S. They have documented the coefficients that are required to evaluate these polynomials for 1000+ species.

The NASA polynomials have the form:

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

R = a1 InT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7

Where R is the universal gas constant, T is the temperature and al, a2, a3, 04, a5, a6, and a7 are the numerical coefficients supplied in the thermodynamic file. The FIRST 7 coefficients
are HIGH-temperature coefficients (above 1000 k) and the SECOND 7 coefficients are LOW temperature coefficients (below 1000 K).

Coding Functions:

1. Specific Heat:
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

Cp= (a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4)*R

Function for specific heat

function Cp = Specificheat (a1, a2, a3, a4, a5, a8, a9, a10, a11, a12, Temp, R, global_medium_temp)
% Comparing with global temperature
% If less than global temperature - minimum coefficients are taken
if Temp <= global_medium_temp
Cp = R*(a8 + (a9*Temp) + (a10*(Temp).^2)+ (a11*(Temp).^3) + (a12*(Temp).^4));
else
Cp = R*(a1 + (a2*Temp) + (a3*(Temp).^2) + (a4* (Temp).^3) + (a5*(Temp).^4));
% Else maximum coefficients are taken
end

2. 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.

S = (a1 lnT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7)*R

function S = Entropy_calc(a1, a2, a3, a4, a5, a7, a8, a9, a10, a11, a12, a14, Temp,global_medium_temp)
% Comparing with global temperature
% If less than global temperature - minimum coefficients are taken
if Temp <= global_medium_temp
S = R.*(a8.*log (Temp) + (a9.*Temp) + (a10*((Temp).^2)/2) + (a11*((Temp).^3)/3) + (a12*( (Temp).^4)/4)+ a14);
else
% Else - maximum coefficients are taken
S = R.*(a1.*log (Temp) + (a2.*Temp)+ (a3*((Temp).^2)/2) + (a4*((Temp).^3)/3) + (a5* ((Temp).^4)/4)+ a7);
end

3. Enthalpy:

Enthalpy, a property of a thermodynamic system, is equal to the system's internal energy plus the product of its pressure and volume. In a system enclosed so as to prevent the mass transfer, for processes at constant pressure, the heat absorbed or released equals the change in enthalpy.

H = (a1 + a2 T/2 + a3 T^2 /3 + a4 T^3 / 4 + a5 T^4 /5 + a6/T)*R*T

% Enthalpy Calculation function
function H = Enthalpy(a1, a2, a3, a4, a5, a6, a8, a9, a10, a11, a12, a13, Temp, R, global_medium_temp)
% Comparing with global temperature
% If less than global temperature - minimum coefficients are taken
if Temp <= global_medium_temp
H = R*Temp.*(a8 + ((a9*Temp)/2) + ((a10*(Temp).^2)/3) + ((a11*(Temp).^3)/4)+((a12*(Temp).^4)/5)+((a13./Temp)));
else
% Else maximum coefficients are taken
H = R*Temp.* (a1 + ((a2*Temp)/2) + ((a3*(Temp).^2) /3) + ((a4*(Temp).^3)/4)+((a5*(Temp).^4/5)+((a6./Temp))));
end

4. Molecular Weight:

The molecular mass is the mass of a given molecule: it is measured in daltons. Different molecules of the same compound may have different molecular masses because they contain different isotopes of an element.

%Molecular Weight Calculation function

function Molecularweight = Molecular_weight(Species_SP)
%Standard - Atomic Name 
atomic_name = ['H','C','O','N','A'];
%Hydrogen,Carbon,Oxygen,Nitrogen,Argon

fprintf('Molecular Weight of');
fprintf(Species_SP);
fprintf(' : ');
% Standard Atomic Number 
Atomic_Weight = [1.008 12.011 15.999 14.007 39.948];
% Starting Range of Molecule
Molecularweight = 0;
% Adding Atomic number to obtain 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 = str2double(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

MAIN CODE

clear all
close all
clc

% Make directory
mkdir 'Desktop NASA File Parsing'
%Reading Data file- NASA Thermodynamic
% Open a file & 'r' Open the file for reading
f1=fopen ('THERMO.dat', 'r');
gt=fgetl(f1);
% Obtain Temperatures - by Reading Header Information
global_temp=fgetl(f1);
% Spliting
t=strsplit(global_temp) ;
% Converting string cells to numbers using stranum,
global_low_temp=str2num(t{2});
global_medium_temp=str2num(t{3});
global_high_temp=str2num(t{4});
% Temperature range & Gas Constant
% Temperature Range assumed
Temp=linspace (global_low_temp, global_high_temp, 1000) ;
% Universal Gas Constant J/ (mol-K)
R = 8.314;
%skipping three comment lines in the file

for i = 1:3
gt=fgetl(f1);
end
for i = 1:53
gt=fgetl(f1);
A=strsplit(gt,' ');
Species_SP=(A{1});
line_1=fgetl(f1);
a=strfind(line_1,'E');
% Extracting higher and lower temperature coefficients
%% Extracting FIRST 7 coefficients i.e HIGH-temperature coefficients
a1=line_1(1:a(1)+3);
a1=str2num(a1);   %1
a2=line_1(a(1)+4:a(2) +3);
a2=str2num(a2);   %2
a3=line_1(a(2)+4:a(3) +3);
a3=str2num(a3);    %3
a4=line_1(a(3)+4:a(4) +3);
a4=str2num(a4);   %4
a5=line_1(a(4)+4:a(5) +3);
a5=str2num(a5);    %5


line_2=fgetl(f1);
b=strfind(line_1,'E');
a6=line_2(1:b(1)+3);
a6=str2num(a6);   %6
a7=line_2(a(1)+4:b(2)+3);
a7=str2num(a7);     %7

%Extracting SECOND 7 coefficinets -LOW-temperature coeficients
a8=line_2 (b(2)+4:b(3)+3);
a8=str2num(a8);    %8
a9=line_2 (b(3)+4: b (4)+3);
a9=str2num(a9);   %9
a10=line_2 (b(4) +4:b (5)+3);
a10=str2num(a10); %10

line_3=fgetl(f1);
c=strfind (line_1,'E');
a11=line_3(1:c(1)+3);
a11=str2num(a11);   %11
a12=line_3(c(1)+4:c(2)+3);
a12=str2num(a12);   %12
a13=line_3(c(2)+4:c(3)+3);
a13=str2num(a13) ;  %13
a14=line_3(c(3)+4:c(4)+3);
a14=str2num(a14);   %14 

% 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
Molecular_list = fopen('Molecular weight.txt','a+');
fprintf(Molecular_list,'Molecular weight of %s and %d n',Species_SP,Molecularweight);
fclose(Molecular_list);

%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 ( 'Varying Specific Heat with Temperature Range of %s',Species_SP));
grid on
savesas(figure (1), 'Specific Heat.jpg');
% Plot 2, Enthalpy (H) vs Temperature Range (Temp)
figure(2)
plot (Temp, H, 'linewidth',2,'color','g');
xlabel( 'Temperature [K]');
ylabel ('Enthalpy in [kJ/ (mol)]') ;
title(sprintf('Varying Enthalpy with Temperature Range of %s',Species_SP));
grid on
saveas(figure (2) ,'Enthalpy.jpg');

% Plot 3, Entropy (S) vs Temperature Range (Temp)
figure(3)
plot (Temp,S, 'linewidth', 2, 'color', 'b');
xlabel ( 'Temperature [K]') ;
label( 'Entropy in (kJ/ (mol-K)]');
title(sprintf( 'Varying Entropy with 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

'clear all' - clear all the variables in the workspace,

'close all' - closes all the plots and figures and

'clc' - clears the command window.
Making a directory using 'mkdir' with location address.
After that, reading the text file 'Thermo.dat' file and opening it by using ‘fopen’ and 'r' as reading command.
For reading the line in a text file using ‘fgetl’.
To obtain the temperature range from the Thermo.dat file, using ‘strsplit’ to split the temperature into a particular string cell.
Used ‘str2num’ for converting string cell to numbers.
Divided the temperature into a string and then converted to a number so that temperature values can be used separately to solve the NASA thermodynamics.
Using ‘for & end’ loops to skip the next three comment lines in the text file.
Reading all the 53 species to plot specific heat(cp), enthalpy(H), and entropy(S), with temperature values, and to find the molecular weights.
Used ‘for’ loop to read 53 species.
Used ‘strfind’ to find string.
Using 'linspace' specified the temperature range from low to high.
Adding universal gas constant value ‘R’.
Extracted higher and lower temperature coefficient by reading line by line from the species.
Calling specific heat, enthalpy, entropy function codes.
Calculating Molecular weight.

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

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

Output:

>> Molecular_weight AR
Molecular Weight of AR : molecular weight AR :3.994800e+01 

ans =

   39.9480

>> Molecular_weight C
Molecular Weight of C : molecular weight C :1.201100e+01 

ans =

   12.0110

>> Molecular_weight C2H
Molecular Weight of C2H : molecular weight C2H :2.503000e+01 

ans =

   25.0300

>> Molecular_weight C2H2
Molecular Weight of C2H2 : molecular weight C2H2 :2.603800e+01 

ans =

   26.0380

>> Molecular_weight C2H3
Molecular Weight of C2H3 : molecular weight C2H3 :2.704600e+01 

ans =

   27.0460

>> Molecular_weight C2H4
Molecular Weight of C2H4 : molecular weight C2H4 :2.805400e+01 

ans =

   28.0540

>> Molecular_weight C2H5
Molecular Weight of C2H5 : molecular weight C2H5 :2.906200e+01 

ans =

   29.0620

>> Molecular_weight C2H6
Molecular Weight of C2H6 : molecular weight C2H6 :3.007000e+01 

ans =

   30.0700

>> Molecular_weight C3H7
Molecular Weight of C3H7 : molecular weight C3H7 :4.308900e+01 

ans =

   43.0890

>> Molecular_weight C3H8
Molecular Weight of C3H8 : molecular weight C3H8 :4.409700e+01 

ans =

   44.0970

>> Molecular_weight CH
Molecular Weight of CH : molecular weight CH :1.301900e+01 

ans =

   13.0190

>> Molecular_weight CH2
Molecular Weight of CH2 : molecular weight CH2 :1.402700e+01 

ans =

   14.0270

>> Molecular_weight CH2(S)
Molecular Weight of CH2(S) : molecular weight CH2(S) :1.402700e+01 

ans =

   14.0270

>> Molecular_weight CH2CHO
Molecular Weight of CH2CHO : molecular weight CH2CHO :4.304500e+01 

ans =

   43.0450

>> Molecular_weight CH2CO
Molecular Weight of CH2CO : molecular weight CH2CO :4.203700e+01 

ans =

   42.0370

>> Molecular_weight CH2O
Molecular Weight of CH2O : molecular weight CH2O :3.002600e+01 

ans =

   30.0260

>> Molecular_weight CH2OH
Molecular Weight of CH2OH : molecular weight CH2OH :3.103400e+01 

ans =

   31.0340

>> Molecular_weight CH3
Molecular Weight of CH3 : molecular weight CH3 :1.503500e+01 

ans =

   15.0350

>> Molecular_weight CH3CHO
Molecular Weight of CH3CHO : molecular weight CH3CHO :4.405300e+01 

ans =

   44.0530

>> Molecular_weight CH3O
Molecular Weight of CH3O : molecular weight CH3O :3.103400e+01 

ans =

   31.0340

>> Molecular_weight CH3OH
Molecular Weight of CH3OH : molecular weight CH3OH :3.204200e+01 

ans =

   32.0420

>> Molecular_weight CH4
Molecular Weight of CH4 : molecular weight CH4 :1.604300e+01 

ans =

   16.0430

>> Molecular_weight CN
Molecular Weight of CN : molecular weight CN :2.601800e+01 

ans =

   26.0180

>> Molecular_weight CO
Molecular Weight of CO : molecular weight CO :2.801000e+01 

ans =

   28.0100

>> Molecular_weight CO2
Molecular Weight of CO2 : molecular weight CO2 :4.400900e+01 

ans =

   44.0090

>> Molecular_weight H
Molecular Weight of H : molecular weight H :1.008000e+00 

ans =

    1.0080

>> Molecular_weight H2
Molecular Weight of H2 : molecular weight H2 :2.016000e+00 

ans =

    2.0160

>> Molecular_weight H2CN
Molecular Weight of H2CN : molecular weight H2CN :2.803400e+01 

ans =

   28.0340

>> Molecular_weight H2O
Molecular Weight of H2O : molecular weight H2O :1.801500e+01 

ans =

   18.0150

>> Molecular_weight H2O2
Molecular Weight of H2O2 : molecular weight H2O2 :3.401400e+01 

ans =

   34.0140

>> Molecular_weight HCCO
Molecular Weight of HCCO : molecular weight HCCO :4.102900e+01 

ans =

   41.0290

>> Molecular_weight HCCOH
Molecular Weight of HCCOH : molecular weight HCCOH :4.203700e+01 

ans =

   42.0370

>> Molecular_weight HCN
Molecular Weight of HCN : molecular weight HCN :2.702600e+01 

ans =

   27.0260

>> Molecular_weight HCNN
Molecular Weight of HCNN : molecular weight HCNN :4.103300e+01 

ans =

   41.0330

>> Molecular_weight HCNO
Molecular Weight of HCNO : molecular weight HCNO :4.302500e+01 

ans =

   43.0250

>> Molecular_weight HCO
Molecular Weight of HCO : molecular weight HCO :2.901800e+01 

ans =

   29.0180

>> Molecular_weight HNCO
Molecular Weight of HNCO : molecular weight HNCO :4.302500e+01 

ans =

   43.0250

>> Molecular_weight HNO
Molecular Weight of HNO : molecular weight HNO :3.101400e+01 

ans =

   31.0140

>> Molecular_weight HO2
Molecular Weight of HO2 : molecular weight HO2 :3.300600e+01 

ans =

   33.0060

>> Molecular_weight HOCN
Molecular Weight of HOCN : molecular weight HOCN :4.302500e+01 

ans =

   43.0250

>> Molecular_weight N
Molecular Weight of N : molecular weight N :1.400700e+01 

ans =

   14.0070

>> Molecular_weight N2
Molecular Weight of N2 : molecular weight N2 :2.801400e+01 

ans =

   28.0140

>> Molecular_weight N2O
Molecular Weight of N2O : molecular weight N2O :4.401300e+01 

ans =

   44.0130

>> Molecular_weight NCO
Molecular Weight of NCO : molecular weight NCO :4.201700e+01 

ans =

   42.0170

>> Molecular_weight NH
Molecular Weight of NH : molecular weight NH :1.501500e+01 

ans =

   15.0150

>> Molecular_weight NH2
Molecular Weight of NH2 : molecular weight NH2 :1.602300e+01 

ans =

   16.0230

>> Molecular_weight NH3
Molecular Weight of NH3 : molecular weight NH3 :1.703100e+01 

ans =

   17.0310

>> Molecular_weight NNH
Molecular Weight of NNH : molecular weight NNH :2.902200e+01 

ans =

   29.0220

>> Molecular_weight NO
Molecular Weight of NO : molecular weight NO :3.000600e+01 

ans =

   30.0060

>> Molecular_weight NO2
Molecular Weight of NO2 : molecular weight NO2 :4.600500e+01 

ans =

   46.0050

>> Molecular_weight O
Molecular Weight of O : molecular weight O :1.599900e+01 

ans =

   15.9990

>> Molecular_weight O2
Molecular Weight of O2 : molecular weight O2 :3.199800e+01 

ans =

   31.9980

>> Molecular_weight OH
Molecular Weight of OH : molecular weight OH :1.700700e+01 

ans =

   17.0070

RESULTS

Specific Heat, Entropy and Enthalpy Plots for O2, N2 and Co2 species :

1. Plots for O2

a. Variation of Specific heat with the temperature of O2

b. Variation of Entropy with the temperature of O2

c. Variation of Enthalpy with the temperature of O2

MOLECULAR WEIGHTS

Conclusion :

Using MATLAB code, the necessary coefficients required to evaluate various thermodynamic parameters for different species are extracted from NASA thermodynamic data file. The specific heat, enthalpy, and entropy for all the species are calculated and plotted against the local temperature range. The molecular weight is calculated for all the species.

References :

https://in.mathworks.com/help/matlab/ref/mkdir.html

https://in.mathworks.com/matlabcentral/answers/317006-create-new-folders-by-mkdir

https://in.mathworks.com/help/matlab/ref/saveas.html

https://in.mathworks.com/matlabcentral/answers/17580-how-to-define-a-path-in-saveas-command

https://in.mathworks.com/help/matlab/referencelist.html?type=function&category=printing-and-exporting

https://in.mathworks.com/help/matlab/ref/msgbox.html

https://in.mathworks.com/matlabcentral/answers/151692-display-value-in-message-box

https://in.mathworks.com/help/matlab/ref/fprintf.html

https://in.mathworks.com/help/matlab/ref/sprintf.html

http://combustion.berkeley.edu/gri-mech/data/nasa_plnm.html

https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19940013151.pdf