Parsing NASA thermodynamic data Using MATLAB

Aim: -

To write a code in Matlab to parse the NASA thermodynamic data file and then calculate Thermodynamic properties, Molecular Weight and Plot nature of the thermodynamic properties of various species.

Objectives: -

1. To write a function to extract the coefficients of the species
2. To write functions that calculate the enthalpy, entropy and specific heat using the NASA polynomials for the given species.
3. To write a function that calculates the molecular weight of each species.
4. To write a function that plots the Cp, Enthalpy and Entropy of an input species.

Introduction: -

Parsing is the process of breaking a data block into smaller chunks by following a set of rules, so that it can be more easily interpreted, managed, or transmitted by a computer. NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H, and S using polynomials. They have also documented the co-efficient that are required to evaluate these polynomials for 1000+ species.

Function to locate the coefficients and store the values in a variables: -

• The input for the below function is the the three strings.
• Next the first string is slected and the location of 'E' is identified in the string.
• After that for the first coefficient, the characters from the start till the position of 'E' and the three characters after the character 'E' are stored in the first coefficient.
• Then for the second a similar process is used but here the starting character for the selection of is the first location of E in the string.
• In the same way all other coefficients are located.
• As the coefficients are in the form of strings, they are converted into numbers by using str2num command.
• In the same way the second and the third strings are used to locate the coefficients.
• After identfying the coefficients they are given as the output of the function.

PROGRAM 1

%Function to identify the coefficient values.
function[ah1,ah2,ah3,ah4,ah5,ah6,ah7,al1,al2,al3,al4,al5,al6,al7] = coefficients(str1,str2,str3)
    %Identifying the location of E in the string
    n = findstr(str1,'E');
    %Storing the values in the left of E and the right of E
    a1 = str1(1:n(1)+3);
    a2 = str1(n(1)+4:n(2)+3);
    a3 = str1(n(2)+4:n(3)+3);
    a4 = str1(n(3)+4:n(4)+3);
    a5 = str1(n(4)+4:n(5)+3);
    %Convertion of the string to a number
    ah1 = str2num(a1);
    ah2 = str2num(a2);
    ah3 = str2num(a3);
    ah4 = str2num(a4);
    ah5 = str2num(a5);
    
    n1 = findstr(str2,'E');
    a6 = str2(1:n1(1)+3);
    a7 = str2(n1(1)+4:n1(2)+3);
    a8 = str2(n1(2)+4:n1(3)+3);
    a9 = str2(n1(3)+4:n1(4)+3);
    a10 = str2(n1(4)+4:n1(5)+3);
    ah6 = str2num(a6);
    ah7 = str2num(a7);
    al1 = str2num(a8);
    al2 = str2num(a9);
    al3 = str2num(a10);
    
    n2 = findstr(str3,'E');
    a1 = str3(1:n2(1)+3);
    a2 = str3(n2(1)+4:n2(2)+3);
    a3 = str3(n2(2)+4:n2(3)+3);
    a4 = str3(n2(3)+4:n2(4)+3); 
    al4 = str2num(a1);
    al5 = str2num(a2);
    al6 = str2num(a3);
    al7 = str2num(a4);
end

Function for calculating molecular weight: -

• First of all the individual elements (Elements) are defined along with there atomic weights (atomic_weight).
• Then the based upon the input elements obtained, the elements are compared with the defined elements by using 'strcmp' command.
• If there is a repetition of an element, ‘str2double’ is used to find that and (n-1) is multiplied to atomic_wt of the element to calculate molecular weight.

• %Function for calculating molecular weight
  function [w] = molecular_weight(element_name)
   
      Elements = {'H' 'C' 'N' 'O' 'AR'};
   
      %atomic weight of each elements
      atomic_weight = [1 12 14 16 40];                                               
      w = 0;                                                                        
   
      for i = 1:length(element_name) 
          for j = 1:length(Elements)  
              %Comparing the input element with the defined elements
              if strcmp(element_name(i), Elements(j)) 
                  w= w + atomic_weight(j);     
   
                  %storing the position of matched element  
                  position = j;                                                           
   
              end
   
          end    
          n = str2double(element_name(i));
              if n > 1                                               
                  w = w + atomic_weight(position)*(n-1); 
              end
     end
  end

  Function to calculate Specific Heat, Enthalpy and Entropy: -

  The input for this function is all the 14 coefficients and all the temperatures.

  The local temperature is checked whether it is less than or equal to the global mid temperature.

  If it is true then the first set is executed by using the first 7 temperature coefficients or else the second set of codes are executed by using the last 7 temperature coefficients.

  The formulas defined are given below:

  %Function to calculate Cp,H and S
  function[Cp,H,S] = thermodynamic_func(ah1,ah2,ah3,ah4,ah5,ah6,ah7,al1,al2,al3,al4,al5,al6,al7,R,T,Global_low_temperature,Global_high_temperature,Global_mid_temperature)
       if T <= Global_mid_temperature
         Cp = R*(ah1 + (ah2*T) + (ah3*T.^2) + (ah4*T.^3) + (ah5*T.^4));
         H = R.*T.*(ah1 + ah2.*(T./2) + ah3.*((T.^2)./3) + ah4.*((T.^3)./4)+ ah5.*((T.^4)./5) + ah6./T);
         S = R*(ah1*log(T) + (ah2*T) + ((ah3*T.^2)/2) + ((ah4*T.^3)/3) + ((ah5*T.^4)/4) + ah7);
         
       else
          
          Cp = R*(al1 + (al2.*T) + (al3.*T.^2) + (al4.*T.^3) + (al5.*T.^4));
          H = R.*T.*(al1 + al2.*(T./2) + al3.*((T.^2)./3) + al4.*((T.^3)./4)+ al5.*((T.^4)./5) + al6./T);
          S = R*(al1.*log(T) + (al2.*T) + ((al3.*T.^2)/2) + ((al4.*T.^3)/3) + ((al5.*T.^4)/4) + al7);
  end

  Function to plot graphs and save them in seperate folders: -

  • In this function the idetifies the current directory by using the command 'pwd'.
  • Then it creates a folder based upon the input variable 'Element_name'.
  • Then it plots Specific heat vs temperature plot.
  • In the program while defining the title a %s is used to print the Element name.
  • In the same way as above all the plots are plotted.
  • Then they are saved in the created folder.
  • Using the ‘cd’ command the folder is changed back to the original folder.
  • The function runs for 53 times through the ‘for’ loop and it creates separate folders for each species and saves all the plots in the folders

  %Function to plot graphs for all elements
  function plotting(Cp,H,S,temp,Element_Name)
      %Identify current directory
      cd(pwd);
      dir_name = Element_Name;
      %Creating a folder based on the Element name 
      mkdir(dir_name);
  
      %Plotting temperature vs specific heat 
      figure(1)
      plot(temp,Cp,'lineWidth',3,'color','r');
      title(sprintf('Specific Heat vs Temperature for %s:',Element_Name));
      xlabel('Temperature');
      ylabel('Specific Heat'); 
  
      %Plotting temperature vs Enthalpy
      figure(2)
      plot(temp,H,'linewidth',3,'color','g');
      title(sprintf('Enthalpy vs Temperature for %s:',Element_Name));
      xlabel('Temperature');
      ylabel('Enthalpy'); 
  
      %Plotting temperature vs Entropy
      figure(3)
      plot(temp,S,'linewidth',3,'color','b');
      title(sprintf('Entropy vs Temperature for %s:',Element_Name));
      xlabel('Temperature');
      ylabel('Entropy'); 
      cd(dir_name);
  
      %Saving the plots
      saveas(1,'Cp_vs_temp.png');
      saveas(2,'H_vs_temp.png');
      saveas(3,'S_vs_temp.png'); 
  
      cd ..;
  
  end

  Writing molecular weights of elements to a file: -

  • The function takes Element name and there molecular weights as input.
  • Then a new file is created with name 'Molecular_Weight.dat' and it is opened in append mode denoted by 'a+' so that we can write the molecular weight of all elements in a single file.
  • Using the ‘fprintf’ command the name of the element along with their molecular weight is written to the file.
  • In the fprintf command %s and %f is used to denote string and number.
  • Finally the file is closed.

  %Writing molecular weight of elements to a file   
  function writing_to_file(Element_Name,weight)   
   
  dir_name = Element_Name;  
  mol_file = fopen('Molecular_Weight.dat','a+');  
  fprintf(mol_file, 'Molecular weight of %s is %f \n', dir_name,weight); 
  fclose(mol_file);  
   
  end

  • Main Program:-
  • The file is opened by using 'fopen' command and then stored in a cell named nasa.
  • Then data is read by using 'textscan' command and stored in x.
  • Then the total no of rows identified by using the command 'length(x)'.
  • After that the file is closed and opened again.
  • After that by using 'fgetl' the first line is read and stored in string First_line.
  • Then line two is read in the same way and a string split(strsplit) command is used to split the string by using the delimeter ' '.
  • Then the temperature values are stored as Global_low, mid and High.
  • After that the local temperature range is defined byusing 'linspace' command.
  • Then R is defined.
  • The comments are ignored by using a for loop.
  • After that the no of times the loop has to executed is calculated.
  • Then the first line is taken split.
  • Then first character is stored in the variable name Element_name.
  • After that the next three lines are given as input to function coefficients.
  • Then the functions to calculate Cp,H,S,molecular weights and to plot the graphs are defined.
  • The loop runs based upon the value of the variable Total_no_of_species.
  • After that the file is closed.

  clear all
  close all
  clc
  
  %Opening the Chemkin NASA thermodynamic data file
  nasa = fopen('THERMO.dat','r');
  x = textscan(nasa,'%s',inf,'Delimiter','\n');
  
  %Identifying the total no of rows
  Total_no_of_rows = length(x{1});
  fclose(nasa);
  nasa = fopen('THERMO.dat','r');
  
  %Reading first line
  First_line = fgetl(nasa);
  
  %Getting the global temperature values
  Line2 = fgetl(nasa);
  s = strsplit(Line2,' ');
  Global_low_temperature = str2num(s{2});
  Global_mid_temperature = str2num(s{3});
  Global_high_temperature = str2num(s{4});
  
  �fining the Local Temperature range
  temp = linspace (Global_low_temperature,Global_high_temperature,500);
  
  %universal gas constant in J/mol-K
  R = 8.314510;
  
  %Ignoring the comments
  for i =1:3
      fgetl(nasa);
  end

  
  
  �lculating the total number of times loop has to be executed
  %As we know that first five lines are does not represent the species
  %content and the lines from the line END is not required I decreased the
  %vaue by four and then divided by four as each species values are
  %represented by four lines
  Total_no_of_species = ((Total_no_of_rows - 10)/4);
  
  for i = 1:Total_no_of_species
      %Reading the line
      L = fgetl(nasa);
      
      %Storing the element name
      e = strsplit(L);
      Element_Name = e{1};
      
      %Reading the three lines which contains the values of the 7 High temperature and 7 low temperature coefficients
      str1 = fgetl(nasa);
      str2 = fgetl(nasa);
      str3 = fgetl(nasa);
  
      %Getting the coefficient values by using the coefficients function
      [ah1,ah2,ah3,ah4,ah5,ah6,ah7,al1,al2,al3,al4,al5,al6,al7] = coefficients(str1,str2,str3);
      
      �lculating specific heat
      [Cp,H,S] = thermodynamic_func(ah1,ah2,ah3,ah4,ah5,ah6,ah7,al1,al2,al3,al4,al5,al6,al7,R,temp,Global_low_temperature,Global_high_temperature,Global_mid_temperature);
      
      %Plotting the CP H S for the Elements
      plotting(Cp,H,S,temp,Element_Name);
      
      �lculating the Molecular Weights
      Molecular_Weight(i) = molecular_weight(Element_Name);
      disp(Molecular_Weight(i))
      Molecular_Weights = strcat(Molecular_Weight(i),Element_Name)
      
      %Function to write molecular weight to a file
      weight = Molecular_Weight(i);
      writing_to_file(Element_Name, weight);
  end
  fclose(nasa);