Project 1 - Parsing NASA thermodynamic data

|| Parsing NASA Thermodynamic Data and Calculating Thermodynamic properties of various gaseous elements Using MATLAB ||

OBJECTIVE :

To parse NASA Thermodynamic Data and to calculate thermodynamic properties of various gaseous elements.

INTRODUCTION :

In 1999, NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Specific Heat, Enthalpy and Entropy. They have also documented the co-efficients 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
    
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 the NASA thermodynamic file.(http://combustion.berkeley.edu/gri-mech/version30/files30/thermo30.dat)

The first 7 numbers starting on the second line of each species entry (five of the second line and first two of the third line) are the seven coefficients (a1 to a7,respectively) for the high temperature range (above 1000K, the upper boundary is specified on the first line of the species entry). The following seven numbers are the coefficients (a1 to a7,respectively) for the low temperature range (below 1000K, the lower boundary is specified on the first line of the species entry).

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

MAIN PROGRAM :

The main program goes like this :

clear all
close all
clc

%=====Reading the 'Thermo.dat' file=====%
file = fopen('Thermo.dat','r');
Name = fgetl(file);
tline = fgetl(file);
A = strsplit(tline,' ');
A(:,1) = [];
global_low_temperature = str2num(A{1});
global_mid_temperature = str2num(A{2});
global_high_temperature = str2num(A{3});

%=====Reading comments=====%
comment_1 = fgetl(file);
comment_2 = fgetl(file);
comment_3 = fgetl(file);

%=====Begining of main loop=====%
for i = 1:53
  species = fgetl(file);
  species_splitted = strsplit(species,' ');
  species_name = species_splitted{1};
  
  local_low_temperature = str2num(species_splitted{end-3});
  local_high_temperature = str2num(species_splitted{end-2});
  local_mid_temperature = str2num(species_splitted{end-1});
  
  
  %==Reading 1st line
  Line_1 = fgetl(file);
  %==Identifying 'E' character so as to differentiate coefficients
  L1 = strfind(Line_1,'E');
  %==Extracting the coefficients based on the position of 'E'
  a1 = str2num(Line_1(1:(L1(1)+3)));
  a2 = str2num(Line_1((L1(1)+4):(L1(2)+3)));
  a3 = str2num(Line_1((L1(2)+4):(L1(3)+3)));
  a4 = str2num(Line_1((L1(3)+4):(L1(4)+3)));
  a5 = str2num(Line_1((L1(4)+4):(L1(5)+3)));
  
  %==Reading the second coefficient line
  Line_2 = fgetl(file);
  L2 = strfind(Line_2,'E');
  a6 = str2num(Line_2(1:(L2(1)+3)));
  a7 = str2num(Line_2((L2(1)+4):(L2(2)+3)));
  a8 = str2num(Line_2((L2(2)+4):(L2(3)+3)));
  a9 = str2num(Line_2((L2(3)+4):(L2(4)+3)));
  a10 = str2num(Line_2((L2(4)+4):(L2(5)+3)));
  
  %==Reading the third coefficient line
  Line_3 = fgetl(file);
  L3 = strfind(Line_3,'E');
  a11 = str2num(Line_3(1:(L3(1)+3)));
  a12 = str2num(Line_3((L3(1)+4):(L3(2)+3)));
  a13 = str2num(Line_3((L3(2)+4):(L3(3)+3)));
  a14 = str2num(Line_3((L3(3)+4):(L3(4)+3)));
  
  %==Function to calculate cp, h and s
  T = linspace(local_low_temperature,local_high_temperature,1000);
  [cp h s] = coeff_calc(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,local_mid_temperature,T);
  
  %==Function to calculate molecular weight
  [mol_wt] = molecular_weight(species_name);
  
  
  %==Creating folders for each speccies
  mkdir(['C:\\Users\\Ayush Umang\\Desktop\\M scripts\\NASA Thermodynamic data\\Species\\',species_name]);
  
  %==Plotting graphs
  figure(1)
  plot(T,s,'linewidth',3,'color','r')
  xlabel('Temperature')
  ylabel('Entropy')
  title(species_name)
  
  figure(2)
  plot(T,h,'linewidth',3,'color','b')
  xlabel('Temperature')
  ylabel('Enthalpy')
  title(species_name)
  
  figure(3)
  plot(T,cp,'linewidth',3,'color','m')
  xlabel('Temperature')
  ylabel('Specific Heat')
  title(species_name)
  
  cd('C:\\Users\\Ayush Umang\\Desktop\\M scripts\\NASA Thermodynamic data\\Species\\')
  cd(species_name)
  saveas(1,'T-s.png')
  saveas(2,'T-h.png')
  saveas(3,'T-cp.png')
  
  cd('C:\\Users\\Ayush Umang\\Desktop\\M scripts\\NASA Thermodynamic data')
  
endfor
fclose(file)

PROGRAM TO CALCULATE ENTHALPY, ENTROPY AND SPECIFIC HEAT :

function[cp h s] = coeff_cal(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,local_mid_temperature,T)

R = 8.314; %==Universal gas constant
for j = 1:length(T)
  if T(j)>local_mid_temperature
    cp(j) = (a1 + a2.*T(j) + a3.*T(j).^2 + a4.*T(j).^3 + a5.*T(j).^4).*R;
    h(j) = (a1 + (a2.*T(j))./2 + (a3.*T(j).^2)./3 + (a4.*T(j).^3)./4 + (a5.*T(j).^4)./5 + a6./T(j)).*R.*T(j);
    s(j) = (a1*log(T(j)) + a2*T(j) + (a3*T(j)^2)/2 + (a4*T(j)^3)/3 + (a5*T(j)^4)/4 + a7)*R;
  else
    cp(j) = (a8 + a9.*T(j) + a10.*T(j).^2 + a11.*T(j).^3 + a12.*T(j).^4).*R;
    h(j) = (a8 + (a9.*T(j))./2 + (a10.*T(j).^2)./3 + (a11.*T(j).^3)./4 + (a12.*T(j).^4)./5 + a13./T(j)).*R.*T(j);
    s(j) = (a8*log(T(j)) + a9*T(j) + (a10*T(j)^2)/2 + (a11*T(j)^3)/3 + (a12*T(j)^4)/4 + a14)*R;
  endif
endfor

endfunction

PROGRAM TO CALCULATE MOLECULAR WEIGHTS :

function [mol_wt] = molecular_weight(target_species)
  elements = {'H','C','N','O','AR'};
  elements_wt = [1,12,14,16,40];
  species = upper(target_species);
  mol_wt = 0;
  
  if (strcmp(species,elements(5)))
    mol_wt = 40;
   else
    for i = 1:length(species)
      for j = 1:length(elements)
        if (strcmp(species(i),elements(j)))
          mol_wt = mol_wt + elements_wt(j);
          v = j;
        endif
      endfor
      n = str2double(species(i));
      if (n>1)
        mol_wt = mol_wt+elements_wt(v)*(n-1);
      endif
    endfor
  endif
  fprintf('molecular weight of %s:- ',target_species)
  fprintf('%d',mol_wt)
  fprintf(' ')
  fprintf('amu ')
  disp(' ')
  
endfunction

EXPLANATION OF THE PROGRAM :

1. The programme starts with 'clc', 'close all' and 'clear all'. 'clc' command is used to clear the command window. 'close all' command closes all the plots and figures which are opened in different windows of MATLAB and 'clear all' command clears all the variables stored in the workplace.
2. After that, the 'Thermo.dat' file provided by NASA, which contains the required data of coefficients, is parsed to the program with the name 'file'.
3. Then the header data is read by using the 'fgetl' command, followed by fetching of global low, mid and high temperature ranges.
4. Comments are read.
5. 'For' loop is used to run the required operation for all the species one by one.
6. Now, the coefficients are fetched and converted to 'num' type data.
7. For the purpose of plotting graphs, a temperature range is set using 'linspace' command with lower limit as 'local_low_temperature', and upper limit as 'local_high_temperature'.
8. After that the function for calculating Enthalpy, Entropy and Specific Heat is called, which is saved in a separate program named 'coeff_calc.m'.
9. Function for calculating the molecular weight is also called.
10. After that, 'mkdir' command is used to make various directories with the species name in the working directory, so that we can save the plots for individual species in thier corresonding directories.
11. Then, three plots are made - Entropy vs Temperature, Enthaply vs Temperature and Specific heat vs Temperature.
12. These plots are saved to their respective directories.
13. Finally, the parsed file is closed using 'fclose' command.

OUTPUTS :

For O2 :

For N2 :

For CO2 :

NOTE : I have only provided plots for O2, N2 and CO2. Outputs for all the species can be accessed here-

https://drive.google.com/open?id=188GP7MeevZROXPnq6M0K3wIXsJpKU-Ah

Output of Molecular Weight function :

REFERENCES :

1. http://combustion.berkeley.edu/gri-mech/data/nasa_plnm.html
2. https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19940013151.pdf

======================Thankyou=====================