BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Project

Parsing NASA thermodynamic data using MATLAB

HHari Prasanth G19 May 2020
MATLAB

AIM
Generate thermodynamic data from numerical coefficients supplied in NASA thermodynamic 'dat' file and plot datas for each species.

OBJECTIVES

  1. Create a function that extracts the 14 co-efficients by reading NASA Polynomial format for CHEMKIN-II
  2. Calculate the enthalpy, entropy and specific heats for all the species in the data file.
  3. Calculate the molecular weight of each species and display it.
  4. Plot the Cp, Enthalpy and Entropy for each species.
  5. Save the plots as images with separate folders for each species

BASIC STUDY

NASA Polynomial format for CHEMKIN-II THERMO

Sample Data in NASA Polynomial format for CHEMKIN-II THERMO

EQUATIONS

The NASA polynomials have the form:

`Cp = (a_1 + a_2*T + a_3*T^2 + a_4*T^3 + a_5*T^4)*R`


`H = R*T*(a_1 + (a_2*T)/2 + (a_3*T^2)/3 + (a_4*T^3)/4 + (a_5*T^4)/5 + a_6/T)`


`S = R*(a_1*log(T) + a_2*T + (a_3*T^2)/2 + (a_4*T^3)/3 + (a_5*T^4)/4 + a_7);`

where `a_1`, `a_2`, `a_3`, `a_4`, `a_5`, `a_6`, and `a_7` are the numerical coefficients supplied in 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

ParsingNASAthermodynamicData

clear all
close all
clc

% change format of integer to long
format long

% Setting folder to use and save data
OutputFolder = ('C:\Users\ADMIN\Desktop\NASA_THERMO2');
ProgramFolder = pwd;

% Open file contains data
fileID = fopen('C:\Users\ADMIN\Desktop\NASA_THERMO\thermo.dat');

% Get Header 
Header = fgetl(fileID);

% Get Temperature ranges
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

% Read Species thermodynamic data

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('********** END reached : Parsing THERMO Data COMPLETED**********')
        break
    end
    fseek(fileID, -(length(next)+1), 0); 
    
    % set folder to OutputFolder
    mkdir(OutputFolder,SpeciesName);
    SpeciesFolder = strcat(OutputFolder,'\',SpeciesName);
    cd(SpeciesFolder);
    
    % 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

clear all
close all

% get Molecular Weight from input
getMolecularWeight()

Function to Extract data from Line 1

% Function to Extract data 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

Function to Extract data from Line 2 and 3

% Function to Extract data from Line 2 and 3

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

Function to Extract data from Line4

% Function to Extract data 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

Calculate thermodynamic properties such as 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

Function to get molecular weight from chemical formula of 'H C N O A' atoms

% Function to get molecular weight from chemical formula of 'H C N O A' atoms

function [molecularWeight] = calcMolecularWeight(SpeciesName)

    % To make search case insenstive
    SpeciesName = upper(SpeciesName);
    
    % Declaring atomic symbols and atomic weight
    AtomicSymbols = {'H','C','N','O','A'};
    atomicWeight = [1 12 14 16 40];
    
    % initializing molecularWeight 0
    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

Function to display molecular weight from input of 'H C N O A' atoms

% Function to get molecular weight from input of 'H C N O A' atoms

function [] = getMolecularWeight()
    while(1)
        % Ask Need to calculate Molecular Weight
        inp = strtrim(input('Need to calculate MolecularWeight of H C N O Ar atoms? (Y / N)','s'));
        inp =upper(inp);
        if(inp=='N')
            disp('********** END **********')
            break
        elseif(inp=='Y')
            inp = strtrim(input('Enter Chemical Formula of H C N O Ar atoms?','s'))
            calcMolecularWeight(inp);
        else
            disp('INVALID INPUT !!!!!!!! Try Again')
        end
    end
end

EXPLANATION

  • STEP 1:
    • Create a function that extracts the 14 co-efficients by reading NASA Polynomial format for CHEMKIN-II
  • STEP 2:
    • Calculate the enthalpy, entropy and specific heats for all the species in the data file.
  • STEP 3:
    • Calculate the molecular weight of each species and display it.
  • STEP 4:
    • Plot the Cp, Enthalpy and Entropy for each species.
  • STEP 5:
    • Save the plots as images with separate folders for each species

RESULTS

MolecularWeight of          O is 16 amu
MolecularWeight of         O2 is 32 amu
MolecularWeight of          H is 1 amu
MolecularWeight of         H2 is 2 amu
MolecularWeight of         OH is 17 amu
MolecularWeight of        H2O is 18 amu
MolecularWeight of        HO2 is 33 amu
MolecularWeight of       H2O2 is 34 amu
MolecularWeight of          C is 12 amu
MolecularWeight of         CH is 13 amu
MolecularWeight of        CH2 is 14 amu
MolecularWeight of     CH2(S) is NaN amu
MolecularWeight of        CH3 is 15 amu
MolecularWeight of        CH4 is 16 amu
MolecularWeight of         CO is 28 amu
MolecularWeight of        CO2 is 44 amu
MolecularWeight of        HCO is 29 amu
MolecularWeight of       CH2O is 30 amu
MolecularWeight of      CH2OH is 31 amu
MolecularWeight of       CH3O is 31 amu
MolecularWeight of      CH3OH is 32 amu
MolecularWeight of        C2H is 25 amu
MolecularWeight of       C2H2 is 26 amu
MolecularWeight of       C2H3 is 27 amu
MolecularWeight of       C2H4 is 28 amu
MolecularWeight of       C2H5 is 29 amu
MolecularWeight of       C2H6 is 30 amu
MolecularWeight of      CH2CO is 42 amu
MolecularWeight of       HCCO is 41 amu
MolecularWeight of      HCCOH is 42 amu
MolecularWeight of       H2CN is 28 amu
MolecularWeight of        HCN is 27 amu
MolecularWeight of        HNO is 31 amu
MolecularWeight of          N is 14 amu
MolecularWeight of        NNH is 29 amu
MolecularWeight of        N2O is 44 amu
MolecularWeight of         NH is 15 amu
MolecularWeight of        NH2 is 16 amu
MolecularWeight of        NH3 is 17 amu
MolecularWeight of         NO is 30 amu
MolecularWeight of        NO2 is 46 amu
MolecularWeight of       HCNO is 43 amu
MolecularWeight of       HOCN is 43 amu
MolecularWeight of       HNCO is 43 amu
MolecularWeight of        NCO is 42 amu
MolecularWeight of         CN is 26 amu
MolecularWeight of       HCNN is 41 amu
MolecularWeight of         N2 is 28 amu
MolecularWeight of         AR is 40 amu
MolecularWeight of       C3H8 is 44 amu
MolecularWeight of       C3H7 is 43 amu
MolecularWeight of     CH3CHO is 44 amu
MolecularWeight of     CH2CHO is 43 amu
********** END reached : Parsing THERMO Data COMPLETED**********
Need to calculate MolecularWeight of H C N O Ar atoms? (Y / N)y
Enter Chemical Formula of H C N O Ar atoms?CH3CH2CH3

inp =

CH3CH2CH3

MolecularWeight of  CH3CH2CH3 is 44 amu
Need to calculate MolecularWeight of H C N O Ar atoms? (Y / N)Y
Enter Chemical Formula of H C N O Ar atoms?NO2

inp =

NO2

MolecularWeight of        NO2 is 46 amu
Need to calculate MolecularWeight of H C N O Ar atoms? (Y / N)n
********** END **********

Plots for o2 species

Plots for N2 species

Plots for Co2 species

CONCLUSION

Thermodynamic data are generated from numerical coefficients supplied in NASA thermodynamic 'dat' file and ploted datas for each species.

ERRORS

Correction

% Function to Extract data 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

REFERENCES

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.