Parsing NASA thermodynamic data Using MATLAB

AIM : -To write a program in MATLAB for parsing NASA thermodynamics data.

OBJECTIVE :

1.      Making of a MATLAB code to extracts 14 coefficients by reading NASA Polynomial format as a .dat file.

2.      From the given DATA file calculating Molecular weight of each species and display it as a result

3.      To plot the Cp, Enthalpy and Entropy for the local temperature range specific for each species save the plots as images with appropriate names and dump them in separate folder.

THEORY : -

Thermodynamic data for individual species are necessary for performing equilibrium chemistry and finite-rate chemical kinetics calculations. These data typically provide the thermodynamic properties (enthalpy, entropy, and specific heat) of a species as a function of temperature. Thermodynamic data for individual species are required in many applications involving chemical reactions. For over 50 years the NASA Lewis (now Glenn) Research Center has been compiling and disseminating thermodynamic data for use in its chemical equilibrium programs. These data, widely used by the thermodynamic community, have grown from 42 species to the current 2000. For calculation purposes, the use of thermodynamic data in the form of simple empirical equations has some obvious advantages. First, it makes tabular interpolations unnecessary; second, it permits analytical integrations; and, third, it condenses all the tabulated information into a few constants that accurately reproduce the many thermodynamic data points.

Another sample :-

EQUATIONS

The NASA polynomials have the form:

Specific heat(Cp): : 

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

Specific heat, ratio of the quantity of heat required to raise the temperature of a body one degree to that required to raise the temperature of an equal mass of water one degree. The term is also used in a narrower sense to mean the amount of heat, in calories, required to raise the temperature of one gram of a substance by one Celsius degree.

Enthalpy(H) :

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

`H = R * T * ( a1 + frac{a2 *T }{2} + (frac{a3 * T^2}{3} +(a4 * T^3)/ 4 + (a5 *T^4)/ 5 + (a6)/T)`
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.

`S = R ⋅ ( a 1 * log ( T ) + a 2 * T + (a3 * T^2)/2 + (a4 *T^3 )/3 + (a5*T^4)/4 + a7 )`

 Here , a1,a2,a3,a4,a5,a6 and a7 are coefficients given in the NASA thermodynamic files.

The first 7 numbers starting on the second line of each species entry are the seven coefficients (a1 through a7, respectively) for the high-temperature range (above 1000 K). The following seven numbers are the coefficients (a1 through a7, respectively) for the low-temperature range (below 1000 K).

 MATLAB PROGRAM 

clear all
close all
clc

% change format of integer to long
format long

% Setting folder to use and save data
OutputFolder = ('Z:\IMPORTANT\local H\GENERAL STUDY\ONLINE LEARNING MODULE\SKILL LYNC\MATLAB FOR MECHANICAL\last week and final submit\PROJECT 1');
ProgramFolder = pwd;

% Open file contains data i.e locaton of parsing file
fileID = fopen('Z:\IMPORTANT\local H\GENERAL STUDY\ONLINE LEARNING MODULE\SKILL LYNC\MATLAB FOR MECHANICAL\last week and final submit\PROJECT 1\THERMO.dat');

%% Get Header (that has to be obselate)
Header = fgetl(fileID);

%Obtaining Temperature ranges according to min,mid and high
TempRange = fgetl(fileID);
lowestT = str2double(TempRange(1,1:10));
commonT = str2double(TempRange(1,11:20));
highestT = str2double(TempRange(1,21:30));
fileComment='';

%% Get data information and comments
while(1)
    comment = fgetl(fileID);
    fileComment = strcat(fileComment,' ',comment);
    next = fgetl(fileID);
    if (next(1,1) ~= '!')
        fseek(fileID, -(length(next)+1), 0);
        break
    end
    fseek(fileID, -(length(next)+1), 0);
end

%% Now extracting data 1 by 1 from the data file

while(1)
    
    % Extract data from line1
    line1 = fgetl(fileID);
    [SpeciesName,Date,AtomicSymbolsAndFormula,Phase,LowTLoc,HighTLoc,CommonTLoc,optAtomicSymbolsAndFormula,Integer] = ExtractLine1(line1);
    
    % Calculate molecular Weight
    molecularWeight = calcMolecularWeight(SpeciesName);
    
    % Extract data from line2
    line2 = fgetl(fileID);
    [a1, a2, a3, a4, a5] = ExtractLine23(line2);
    
    % Extract data from line3
    line3 = fgetl(fileID);
    [a6, a7, a8, a9, a10] = ExtractLine23(line3);
    
    % Extract data from line4
    line4 = fgetl(fileID);
    [a11, a12, a13, a14] = ExtractLine4(line4);
    
    % Inputs to Generate thermodynamic data
    TemperatureRangeLoc = [LowTLoc; HighTLoc; CommonTLoc];
    coefficients = [a1; a2; a3; a4; a5; a6; a7; a8; a9; a10; a11; a12; a13; a14];
    
    % Generate thermodynamic data
    [Cp, H, S, Temperature] = CalcThermoProp(coefficients,TemperatureRangeLoc);
    
    % Check for 'END' to complete parsing
    next = fgetl(fileID);
    if (strcmp(next,'END'))
        disp('***Hurray SUCCESSFULL  : Parsing THERMO Data COMPLETED***')
        break
    end
    fseek(fileID, -(length(next)+1), 0); 
    
    % set folder to OutputFolder
    mkdir(OutputFolder,SpeciesName);
    SpeciesFolder = strcat(OutputFolder,'',SpeciesName);
    cd(SpeciesName);
    
    % plot Temperature vs Cp
    figure(1)
    plot(Temperature,Cp);
    xlabel('Temperature in [K]');
    ylabel('Cp in [J/(mol K)]');
    title(strcat(SpeciesName,'-','Temperature vs Cp'));
    saveas(gcf,'Temperature vs Cp','jpeg')
    
    % plot Temperature vs Enthalpy
    plot(Temperature,H);
    xlabel('Temperature in [K]');
    ylabel('Enthalpy  in [J]');
    title(strcat(SpeciesName,'-','Temperature vs Enthalpy'));
    saveas(gcf,'Temperature vs Enthalpy','jpeg')
    
    % plot Temperature vs Entropy 
    plot(Temperature,S);
    xlabel('Temperature in [K]');
    ylabel('Entropy  in [J/K]');
    title(strcat(SpeciesName,'-','Temperature vs Entropy '));
    saveas(gcf,'Temperature vs Entropy ','jpeg')
    
    % set folder to ProgramFolder
    cd(ProgramFolder);
end

SOME SUB-FUNCTION THAT IS INVOLVED IN THE MAIN FUNCTION :-

• To get data from Line 1.

  %% data extration from Line 1
  
  function [SpeciesName,Date,AtomicSymbolsAndFormula,Phase,LowT,HighT,CommonT,optAtomicSymbolsAndFormula,LineInteger] = ExtractLine1(line)
  	SpeciesName = strtrim(line(1,1:18));
      Date = strtrim(line(1,19:24));
  	AtomicSymbolsAndFormula = strtrim(line(1,25:44));
  	Phase = strtrim(line(1,45));
  	LowT = str2double(line(1,46:55));
  	HighT = str2double(line(1,56:65));
  	CommonT = str2double(line(1,66:75));
  	optAtomicSymbolsAndFormula = strtrim(line(1,76:79));
  	LineInteger = str2double(line(1,80));
  end​

• To get data from Line 2 and Line 3

  %% data extraction from Line2 and Line3
  
  function [a1, a2, a3, a4, a5] = ExtractLine23(line)
      a1 =str2double(line(1,1:15));
      a2 =str2double(line(1,16:30));
      a3 =str2double(line(1,31:45));
      a4 =str2double(line(1,46:60));
      a5 =str2double(line(1,61:75));
  end​

• To get data from Line 4

  %% data extraction from Line4
  
  function [a1, a2, a3, a4] = ExtractLine4(line)
      a1 = str2double(line(1,1:15));
      a2 = str2double(line(1,16:30));
      a3 = str2double(line(1,31:45));
      a4 = str2double(line(1,46:60));
  end​

• finally exacting thermodynamic properties like Cp,H,and S.

  %% function to calculate thermodynamic properties such as Cp, H and S
  
  function [Cp, H, S, Temperature] = CalcThermoProp(coefficients,TemperatureRangeLoc);
  
      % inputs 
      R = 1.987; % The gas constant R in J / mol·K
      
      % set Temperature Range
      LowTLoc = TemperatureRangeLoc(1,1);
      HighTLoc = TemperatureRangeLoc(2,1);
      loopCount = 1;
  
      for i=LowTLoc:10:HighTLoc
          T = i;
          % 
          if (T>1000)
              % high-temperature range (above 1000 K). 
                  a1 = coefficients(1,1);
                  a2 = coefficients(2,1);
                  a3 = coefficients(3,1);
                  a4 = coefficients(4,1);
                  a5 = coefficients(5,1);
                  a6 = coefficients(6,1);
                  a7 = coefficients(7,1);
          else
              % low-temperature range (below 1000 K).
                  a1 = coefficients(8,1);
                  a2 = coefficients(9,1);
                  a3 = coefficients(10,1);
                  a4 = coefficients(11,1);
                  a5 = coefficients(12,1);
                  a6 = coefficients(13,1);
                  a7 = coefficients(14,1);     
          end
          
          % Calculate the enthalpy, entropy and specific heat
          Cp(loopCount) = (a1 + a2*T + a3*T^2 + a4*T^3 + a5*T^4)*R;
          H(loopCount) = R*T*(a1 + (a2*T)/2 + (a3*T^2)/3 + (a4*T^3)/4 + (a5*T^4)/5 + a6/T);
          S(loopCount) = R*(a1*log(T) + a2*T + (a3*T^2)/2 + (a4*T^3)/3 + (a5*T^4)/4 + a7);
           
          Temperature(loopCount) = T;
          loopCount = loopCount + 1;
      end
      
  end​

• To get molecular weight from chemical formula of 'H C N O A' atoms.

  %% Code to Extract molecular weight from chemical formula of 'H C N O A' atoms from DATA file
  
  function [molecularWeight] = calcMolecularWeight(SpeciesName)
  
      % To make search cases depending upon the species
      SpeciesName = upper(SpeciesName);
      
      % Atomic symbols and atomic weight Respectiely
      AtomicSymbols = {'H','C','N','O','A'};
      atomicWeight = [1 12 14 16 40];
      
      % initializing molecularWeight with zero;
      molecularWeight = 0;
      
      % Calculating molecularWeight
      for i = 1:length(SpeciesName)
          
          % Calculating molecularWeight from atom character in chemical formula
          for j = 1:length(AtomicSymbols)
              if strcmp(SpeciesName(i),AtomicSymbols(j))
                  molecularWeight = molecularWeight + atomicWeight(j);
                  x = j; 
                  % Storing last calculated atom character in chemical formula
              end
          end
          
          n = str2double(SpeciesName(i));
          
          % Calculating molecularWeight from atom count in chemical formula
          if (isletter(SpeciesName(i))==0)
              molecularWeight = molecularWeight + atomicWeight(x)*(n-1);
          end
      end
      
      % display Molecular Weight
      disp(sprintf('MolecularWeight of %10s is %d amu',SpeciesName,molecularWeight));
  end​

Arrangement Steps- 

1. To peruse the given information document, fopen order was utilized and fgetl order was utilized to peruse the resulting lines, which were headers, remarks and numerical passages. 
2. To remove the numerical estimations of worldwide and neighborhood temperatures which were put away as string by utilizing fgetl, strsplit order was utilized to part the given line as indicated by spaces gave in the middle of them and were changed over from string to twofold datatype utilizing str2double order and were put away in separate factors. 
3. The absolute number of species present in the given information record was 53, so a for circle was started from 1 to multiple times for perusing every specie and 14 coefficients were available under every specie. 
4. These coefficients were again separated and changed over from string to twofold datatype. 
5. Two separate capacities for atomic weight and explicit warmth, enthalpy, entropy were made which stores the necessary information and equations for figuring. 
6. The capacities made were brought in primary program for additional computations and particular charts were plotted for every specie and put away in singular organizers utilizing compact disc and mkdir order. 
7. At the point when the program is executed sub-atomic load for every specie were shown in the order window and plots of explicit warmth, enthalpy, entropy, w.r.t time for every specie were acquired.

Result and Conclusions :-

Some sample results are as follows : -

For NH3 :-

FOR C : - 

Finally we can conclude that the thermodynamic data was sucessfully extracted and plot of all thermodynamic properties such as Enthalpy,Entropy, and Specific Heat as shown above.

Error :- 

• The following error occured, i.e., the masses of compounds were coming wrong and was due to not storing the output of 'j' into a separate variable, It was corrected by storing j as pos in line 15.

=R⋅T⋅(a1+a2⋅T2+a3⋅T23+a4⋅T34+a5⋅T45+a6T)

Cp=(a1+a2⋅T+a3⋅T2+a4⋅T3+a5⋅T4)⋅R