top of page
Search

Parsing NASA Thermodynamics Data File to Calculate Various Terms Using MATLAB

  • Writer: Mohit Sachdeva
    Mohit Sachdeva
  • Jul 14, 2020
  • 5 min read

Aim:- Parsing NASA Thermodynamics Data File to Calculate Various Terms Using MATLAB.


Software Used:- MATLAB


Objective:-

  • To read “NASA Thermodynamic Data file” and to write a function that extracts the 14 coefficients and calculate the enthalpy, entropy and specific heat for all the species.

  • Calculate the molecular weight of each species and display it in the command window.

  • Plot the Cp, Enthalpy and Entropy for the local temperature range.

  • Save the plots as images with appropriate names and dump them in separate folders for each species.


Introduction:-

Parsing a file means reading in a data stream of some sort and building an in memory model of the semantic content of that data. This aims at facilitating performing some kind of transformation on the data.

For this project we are going to parsing NASA Thermodynamic format for CHEMKIN-II file which contain temperature ranges with 53 species. Thermo.dat file consist of sections containing element data, species data, thermo data, and reaction rate data.



ree


Governing Equation:-


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`

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


here ,

  • R is the universal gas constant (R = 8.314 J/mol. K)

  • T is the temperature (Local Temperatures)

  • Cp is specific heat at constant pressure

  • 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 upper boundary is specified on the first line of the species entry).

  • The next seven numbers are the coefficients for the low-temperature range below 1000 K(, the lower boundary is specified on the first line of the species entry).

  • 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).



Work Flow:-

  1. First we will open and read the file named THERMO.dat.

  2. Then we will extract the header line and global temperature range.

  3. In the given file, lines starting with exclamation signs are just comments so we will skip them in our program using a for loop.

  4. After that, we will extract local temperature range.

  5. Finally comes another for loop with 53 iterations because the numbers of species in the given file are also 53 and all further operations will be done in this loop only with respect to each species.

  6. In this loop, name of species is extracted and the further 3 lines are read to get the position of ‘E’ so that we can extract all 14 coefficients of each species.

  7. First 7 coefficients are the high temperature coefficients while the other 7 are the lowest temperature coefficients.

  8. ‘R’ is the universal gas constant which will be used for further calculations in our made functions.

  9. To calculate enthalpy, entropy, molecular weight and specific heat; we have made separate functions so that our main program does not get complicated

  10. Similarly, we have to plot graphs for each species and save them all in relevant folders as per the species’ names. So for that also we have made a function.

  11. These all functions are then called in the for loop to calculate for each species.

  12. Additionally, I have used some commands to make a txt file to save a list of molecular weights of each species which will be also shown in the command window.

  13. All these files and folders are made automatically using commands, hence reducing the manual effort and time.

NOTE:- A function should be saved as the name of its own because that name will be used to call it in the main program.



Commands used in Program Code:-


fopen(‘filename’,r) - Open and read a file

fgetl - Read line from file

strsplit - Split string or character vector at specified delimiter

Fopen(‘filename’,wt) - open and write in a file

str2num - Convert 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/folder

cd - Change working/current directory

strcmpi - compare string case insensitive

Program Code:-

Specific Heat Function:-

% Function to calculate specific heat(CP)

function CP = Specificheat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,local_med)

% Comparing with local temperature

% If temperature range is greater than local_med, it will take maximum coefficients

if T > local_med

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

% If temperature range is less than local_med, it will take minimum coefficients

else

CP =(b1+b2.*T+b3.*T.^2+b4.*T.^3+b5.*T.^4).*R;

end



Entropy Function:-

% Function to calculate entropy(S)

function S = entropy(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T,R,local_med)

% Comparing with local temperature

% If temperature range is greater than local_med, it will take maximum coefficients

if T > local_med

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

else

% If temperature range is less than local_med, it will take minimum coefficients

S =(b1.*log(T)+b2.*T+b3.*(T.^2)./2+b4.*(T.^3)./3+b5.*(T.^4)./4+b7).*R;

end



Enthalpy Function:-

% Function to calculate enthalpy(H)

function H = Enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T,R,local_med)

% Comparing with local temperature

% If temperature range is greater than local_med, it will take maximum coefficients

if T > local_med

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

else

% If temperature range is less than local_med, it will take minimum coefficients

H =(b1.*T+(b2.*T.^2)./2+(b3.*T.^3)./3+(b4.*T.^4)./4+(b5.*T.^5)./5+b6).*R;

end



Molecular Weight Function:-

% Function to calculate molecular weight(MW)

function MW = Mol_wt(species)

elements=['H','O','C','N','A']; % Standard Atomic Name

atomic_weight=[1.007,15.999,12.010,14.006,39.948]; % Standard Atomic Number

MW=0; % Starting Range of Molecule

% Adding Atomic number to get particular value

for i=1:length(species)

for j=1:length(elements)

%comparing strings case insensitive

if strcmpi(species(i),elements(j))==1

MW=MW+atomic_weight(j);

position=j; %storing position of matched species

end

end

n = str2num(species(i));

if n>1

MW=MW + (atomic_weight(position)*(n-1));

end

end

disp(['Molecular Weight of ',species,' is '])

disp(MW)

end


Plotting Function:-

% Plotting for each species

function Plotting(CP,H,S,T,species)

% Plot 1, Specific heat (CP) vs Temperature Range (T)

figure(1)

plot(T,CP,'linewidth',2,'color','r');

xlabel('Temperature [K]');

ylabel('Specific Heat [kJ]');

title(sprintf('Specific Heat vs Temperature Range of %s',species));

grid on

saveas(figure(1),'Specific Heat.jpg');

% Plot 2, Enthalpy (H) vs Temperature Range (T)

figure(2)

plot(T,H,'linewidth',2,'color','b');

xlabel('Temperature [K]');

ylabel('Enthalpy in [kJ/(mol)]');

title(sprintf('Enthalpy vs Temperature Range of %s',species));

grid on

saveas(figure(2),'Ethalpy.jpg');

% Plot 3, Entropy (S) vs Temperature Range (T)

figure(3)

plot(T,S,'linewidth',2,'color','g');

xlabel('Temperature [K]');

ylabel('Entropy in [kJ/(mol-K)]');

title(sprintf('Entropy vs Temperature Range of %s',species));

grid on

saveas(figure(3),'Entropy.jpg');

end



Main Program Code:-

%program code for parsing nasa thermo file and performing varous

%calculations

close all

clear all % clears workspace

clc % clears command window

%open and read data file

f1=fopen('THERMO.dat','r');

%reading data file first line

header=fgetl(f1)

%global temperature

l1= fgetl(f1); %read second line containing global temperatures

Temp=strsplit(l1,' '); %splitting them after each gap

global_min=str2num(Temp{2});

global_mid=str2num(Temp{3});

global_max=str2num(Temp{4});

%not considering the next 3 lines as they are comments in file

for i=1:3

l2=fgetl(f1);

end

list = fopen('molecular_wt.txt','wt'); %creating txt file to make list of mol wt

for j=1:53 %number of iterations equal to species

%reading next line to get temperatures

l3=fgetl(f1); %Reading line1 of each species

at=strfind(l3,'.'); %to find position of dot

t=strsplit(l3(at-3:end-1),' '); %extracting temperature from line

%extracting local temperature range

local_min=str2num(t{1});

local_med=str2num(t{2});

local_max=str2num(t{3});

T=linspace(local_min,local_max,200);

species=strtok(l3,' '); %extracting species name

%extracting co-efficient for each species containing 3 lines

a=fgetl(f1);

b=fgetl(f1);

c=fgetl(f1);

%finding position of E for each species containing 3 lines

d=strfind(a,'E');

e=strfind(b,'E');

f=strfind(c,'E');

%extracting high temperature coefficients

a1=str2num(a(1:d(1)+3));

a2=str2num(a(d(1)+4:d(2)+3));

a3=str2num(a(d(2)+4:d(3)+3));

a4=str2num(a(d(3)+4:d(4)+3));

a5=str2num(a(d(4)+4:d(5)+3));

a6=str2num(b(1:e(1)+3));

a7=str2num(b(e(1)+4:e(2)+3));

%extracting low temperature coefficients

b1=str2num(b(e(2)+4:e(3)+3));

b2=str2num(b(e(3)+4:e(4)+3));

b3=str2num(b(e(4)+4:e(5)+3));

b4=str2num(c(1:f(1)+3));

b5=str2num(c(f(1)+4:f(2)+3));

b6=str2num(c(f(2)+4:f(3)+3));

b7=str2num(c(f(3)+4:f(4)+3));

%universal gas constant (J/mol.K)

R = 8.3145;

%specific heat,enthalpy,entropy,molecular weight

[CP]=Specificheat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,local_med);

[H] = Enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T,R,local_med);

[S] = entropy(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T,R,local_med);

[MW] = Mol_wt(species);

%creating different folders for respective species

mkdir(species); %make new folder named species

current_dir=pwd; %identify current folder

cd(species) %change current folder

Plotting(CP,H,S,T,species); %creating plots and saving in each folder of species

cd(current_dir); %change current folder

fprintf(list,'Molecular Weight of %s is %f \n\n',species,MW); %writing data in txt file

end

fclose(list); %closing txt file


Output/Result:-

Workspace:-

ree
ree


Command Window:-

ree


Plots/Graphs:-

(O2) Entropy vs Temperature Range

ree

(O2) Enthalpy vs Temperature Range

ree

(O2) Specific Heat vs Temperature Range

ree

(CO2) Entropy vs Temperature Range

ree

(CO2) Enthalpy vs Temperature Range

ree

(CO2) Specific Heat vs Temperature Range

ree

(N2) Entropy vs Temperature Range

ree

(N2) Enthalpy vs Temperature Range

ree

(N2) Specific Heat vs Temperature Range

ree

Screenshot of Saved Folders:-


ree


Conclusion:-

Hence, we have done file parsing and used various functions to solve the calculations.

The use of functions helps to make our main program look neat also.

Here is the link to the folders of each species generated automatically in the program containing plots and also a txt file containing list of molecular weight of each species.


Reference:-




Comments

bottom of page