Project - Parsing NASA thermodynamic data (MATLAB)

Aim: A MATLAB Code to parse the NASA thermodynamic data file and then calculate the thermodynamic properties of various gas species.

NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H, and S. They have also documented the coefficients that are required to evaluate these polynomials for 1000+ species.

1. Write a function that extracts the 14 co-efficient and calculates the enthalpy, entropy, and specific heats for all the species in the data file. (Thermo.dat)
2. Calculating the molecular weight of each species and showing in the command window.
3. Plotting the Cp, Enthalpy, and Entropy for the local temperature range (low temperature: high temperature) specific for each species.
4. Plots for O₂, N_2, and CO₂

The formulae for calculating the NASA Thermodynamic are given

               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 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 next seven numbers are the coefficients for the low-temperature range below 1000 K.
• 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). 

• R is the universal gas constant = 314 J/(mol-K)
• T is the Temperature Range

Thus, for each species, there are 14 coefficients in all. In addition to the polynomial coefficients for each species, the Thermodynamic Data (Thermo.dat txt file) provides the species name, its elemental makeup, and the temperature ranges over which the fits are valid.

The important input to the fitting procedure is specific heat, enthalpy, and entropy as a function of temperature. In addition, the extent of the fit two temperature ranges, the temperature ranges have to be specified. Thus, the common temperature connecting the two ranges is 1000 K,

File Parsing: - File Parsing is a technique to read the information from a particular text file and that information will go to be used by MATLAB. If we know file parsing then we can read data from any other external file and add that data for our calculation. By doing we can perform data transformation. In other words, we can split large data into small data, and we can able to store them so that we can manipulate the data.

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

Important Syntax Used during Coding

Steps Involve in Coding: -

• First of all, before coding the final NASA Thermodynamics, we have to create different functions related to Molecular weight, Read and writing molecular weight, Specific heat, Enthalpy, Entropy & Plotting.
• The function should be created by using if, else, and for loops.
• After creating the function code save the code with an appropriate name so that while calling the functions, we will not get any errors.

Detail Coding of functions

• Specific heat Calculation function code:

The specific heat is the amount of heat per unit mass required to raise the temperature by one degree Celsius.

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

% Calculation of specific heat
 
function CP = Specificheat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp)
         j = 1:length(Temp)
             if Temp (j) <3000
             end
             
% Comparing with global temperature
% If less than global it will take minimum coefficients
 
if  Temp > global_medium_temp
CP =(a1+a2.*Temp+a3.*Temp.^2+a4.*Temp.^3+a5.*Temp.^4).*R;
 
% If temperature is more it will take maximum coefficients
 
else
CP =(a8+a9.*Temp+a10.*Temp.^2+a11.*Temp.^3+a12.*Temp.^4).*R;
 
end

• Enthalpy function code:

A thermodynamic quantity is equivalent to the total heat content of a system. It is equal to the internal energy of the system plus the product of pressure and volume. Enthalpy is similar to energy, but not the same.

H=(a8+(a9.*t./2)+(a10.*(t.^2)./3)+(a11.*(t.^3)./4)+(a12.*(t.^4)./5)+ (a13./t)).*R.*t

% Calculation of enthalpy
 
function H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,Temp,R,global_medium_temp)
            j = 1:length(Temp)
             if Temp (j) <3000
             end
             
% Comparing with global temperature
% If less than global it will take minimum coefficients
 
       if Temp > global_medium_temp
        H =(a1.*Temp+(a2.*Temp.^2)./2+(a3.*Temp.^3)./3+(a4.*Temp.^4)./4+(a5.*Temp.^5)./5+a6).*R;
       else
% If temperature is more it will take maximum coefficients 
           
        H =(a8.*Temp+(a9.*Temp.^2)./2+(a10.*Temp.^3)./3+(a11.*Temp.^4)./4+(a12.*Temp.^5)./5+a13).*R;
end

• Entropy function code:

A thermodynamic quantity representing the unavailability of a system's thermal energy for conversion into mechanical work often interpreted as the degree of disorder or randomness in the system. Entropy is also a measure of the number of possible arrangements the atoms in a system can have.

s=(a8*log(t)+a9.*t+(a10.*(t.^2)./2)+(a11.*(t.^3)./3)+(a12.*(t.^4)./4)+a14).*R

% Calculation of entropy
 
function S = Entropy_calcEntropy_calc(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,Temp,R,global_medium_temp)
         j = 1:length(Temp)
             if Temp (j) <3000
             end
             
% Comparing with global temperature
% If less than global it will take minimum coefficients
 
    if Temp > global_medium_temp
        S =(a1.*log(Temp)+a2.*Temp+a3.*(Temp.^2)./2+a4.*(Temp.^3)./3+a5.*(Temp.^4)./4+a7).*R;
    else
 
% If temperature is more it will take maximum coefficients
        S =(a8.*log(Temp)+a9.*Temp+a10.*(Temp.^2)./2+a11.*(Temp.^3)./3+a12.*(Temp.^4)./4+a14).*R;
end

• Molecular Weight function code: -

Molecular weight is a measure of the sum of the atomic weight values of the atoms in a molecule. Molecular weight is used in chemistry to determine stoichiometry in chemical reactions and equations.

• Writing function for the molecular weight we have to specify a specific atomic number.
• The atomic number should match the atomic name.
• For this calculation here I am using
  • Hydrogen =          1.008
  • Carbon =          12.011
  • Oxygen =          15.999
  • Nitrogen =          14.007
  • Argon =          39.948
• Now provide an initial starting range of molecular weight by 0.
• By using for loop add atomic number with each species to get particular molecular weight.

% Function to calculate molecular weight
 
function Molecularweight = Molecular_weight(Species_SP)
         
         atomic_name = ['H','C','O','N','A'];             % Standard Atomic Name
         % Hydrogen, Carbon, Oxygen, Nitrogen, Argon
         
                fprintf('Molecular Weight of ');
                fprintf(Species_SP);
                fprintf(' : ');
                Atomic_Weight = [1.008 12.011 15.999 14.007 39.948];     % Standard Atomic Number
                    
         Molecularweight = 0;                          % Starting Range of Molecule
         
        % Adding Atomic number to get 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 = str2num(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

• Rear & Write molecular weight function

• This code is basically to create a text file for molecular weight.

% Read & Write Molecular text file
 
function Read_write (Species_SP,Molecularweight)
 
         Molecular_list = fopen('Molecular weight.txt','a+');
         fprintf( Molecular_list,'Molecular weight of %s is %d \n',Species_SP,Molecularweight);
         fclose(Molecular_list);
         
end

• Plotting function code
• In plotting code, we will plot an individual CP, H, S Vs Temperature.

• Also, we will save all the graphs to separate folders.

% Plotting for for each species
 
function 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('Specific Heat vs Temperature Range of %s',Species_SP));
        grid on
        saveas(figure(1),'Specific Heat.jpg');
 
        % Plot 2, Enthalpy (H) vs Temperature Range (Temp)
 
        figure(2)
        plot(Temp,H,'linewidth',2,'color','b');
        xlabel('Temperature [K]');
        ylabel('Enthalpy in [kJ/(mol)]');
        title(sprintf('Enthalpy vs Temperature Range of %s',Species_SP));
        grid on
        saveas(figure(2),'Ethalpy.jpg');
 
        % Plot 3, Entropy (S) vs Temperature Range (Temp)
 
        figure(3)
        plot(Temp,S,'linewidth',2,'color','g');
        xlabel('Temperature [K]');
        ylabel('Entropy in [kJ/(mol-K)]');
        title(sprintf('Entropy vs 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

• Steps Involved in Main Code

• Start writing the program by clear and clean the workspace and command also close the background dialogue box.
• Next Step is to provide default directory by using ‘cd’ to start the program and save all the thing into that directory, here in the program I have used my default location.
• After that, we have to read the text file i.e Thermo.dat file and open it by using ‘fopen’ and r’’.
• Now for reading the line use ‘fgetl’.

                                     Fig 1.1 Detail of Thermo.dat file

• Now next step is to obtain the temperature range from the Thermo.dat file.
• Now by using ‘strsplit’ we have to split the temperature into a particular string cell.
• After that, we have to use ‘str2num’ for converting string cell to numbers.
• Here we divided the temperature into a string and then converted to a number so that we can use separately to solve the NASA thermodynamics.
• Now in the text file, we have to skip the next three lines which we are not useful for coding. So, for that, we are using ‘for & end’ loops.
• After that, we have to read all the 53 species so that we can plot CP,H, and S. Also we can find the molecular weights.
• For reading all those 53 species here I am using ‘for’ loop.
• Now read the first line, here I use O_RL = fgetl.
• For reading specific species use ‘strtok’. This syntax will return to specific species after every reading.
• By ‘strfind’ we will use to find string find from another.
• Now we have to code for reading the local temperature range.
• By using linspace specify the temperature range from low to high and by using 1000 K.
• Also, add universal gas constant value ‘R’.
• Now next step is to extract higher and lower temperature coefficient by reading line by line from the species line.
• Now specify the specific heat, enthalpy, entropy for calling the function code.
• For calculating Molecular weight specify the code.

clear all; close all; clc
 
% Setting default folder on every startup of Matlab 
 
cd 'D:\3D Design & Simulation\Skill lync\Challenges\Parsing NASA thermodynamic data\NASA Thermodynamic\Species'
 
...
 
pwd
 
%  Reading Data file - NASA Thermodynamic
 
f1 = fopen ('THERMO.dat','r');              % Open a file & 'r' Open the file for reading  
gt = fgetl(f1);
 
% Obtaing Temperatures (Reading Header Information)
 
global_temp = fgetl (f1);
temp_range = ('300.000  1000.000  5000.000');                % Spliting
 
A = strsplit(temp_range,' ');
 
% Now convert string cells tonumbers using str2num,
 
global_low_temp = str2num (A{1});                       % Global Low
global_medium_temp = str2num (A{2});                    % Global Medium
global_high_temp = str2num (A{3});                      % Global High
 
% skipping the next three comments in the file
 
for i=3:5
comments = fgetl(f1);
end
 
for i = 1:53
    
    O_RL = fgetl(f1);                               % RL = Read line
    Species_SP = strtok(O_RL);                            % SP = Split
    L_T = strfind (O_RL,'G');
    
    % reading local temperatue range
    
    Avg_temp = O_RL(L_T+1:end-1);
    Avg_temp = str2num(Avg_temp);
    
    % Range of temperature & Gas Constant
 
    Temp = linspace(global_low_temp,global_high_temp,1000);   % Temperature Range assumed
    R = 8.314;                                                % Universal Gas Constant J/(mol-K)
    
    % Extracting higher and lower temperature coeffs
    
    line_1 = fgetl(f1);
    line_2 = fgetl(f1);
    line_3 = fgetl(f1);
    
    A = strfind(line_1,'E');
    B = strfind(line_2,'E');
    C = strfind(line_3,'E');
  
    % Extracting HIGH-temperature coefficients - FIRST 7 coefficients
    
    a1 = str2num(line_1(1:A(1)+3));
    a2 = str2num(line_1(A(1)+4:A(2)+3));
    a3 = str2num(line_1(A(2)+4:A(3)+3));
    a4 = str2num(line_1(A(3)+4:A(4)+3));
    a5 = str2num(line_1(A(4)+4:A(5)+3));
    a6 = str2num(line_2(1:B(1)+3));
    a7 = str2num(line_2(B(1)+4:B(2)+3));
    
    % Extracting LOW-temperature coefficients - SECOND 7 coefficients
    
    a8 = str2num(line_2(B(2)+4:B(3)+3));
    a9 = str2num(line_2(B(3)+4:B(4)+3));
    a10 = str2num(line_2(B(4)+4:B(5)+3));
    a11 = str2num(line_3(1:C(1)+3));
    a12 = str2num(line_3(C(1)+4:C(2)+3));
    a13 = str2num(line_3(C(2)+4:C(3)+3));
    a14 = str2num(line_3(C(3)+4:C(4)+3));
    
    % 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
    Read_write(Species_SP,Molecularweight)
    
    % Creating individual folder for each species and plotting
    Plotting(CP,H,S,Temp,Species_SP)
    
end

Output Results:

Plotting of O2 Graph

                                Fig 1.2 Specific Heat Vs Temperature Range

                                Fig 1.3 Enthalpy Vs Temperature Range

                                       Fig 1.4 Entropy Vs Temperature Range

Plotting of N2 Graph

                                    Fig 1.5 Specific Heat Vs Temperature Range

                                     Fig 1.6 Enthalpy Vs Temperature Range

                                      Fig 1.7 Entropy Vs Temperature Range

Plotting of Co2 Graph

                                     Fig 1.8 Specific Heat Vs Temperature Range

                                  Fig 1.9 Enthalpy Vs Temperature Range

                                    Fig 1.10 Entropy Vs Temperature Range

• Screenshot of the folder location

• Molecular Weight Screenshot