Parsing NASA Thermodynamic Data

Introduction:

Nasa came up with polynomials that can be used to evaluate thermodynamic properties such as CP(Specific Heat), H(Enthalpy), S(Entropy). They have also documented the co_efficients that are required to evaluate these polynomials for 1000+ species.

Aim:

1.To write a program that extracts the 14 coefficients of 53 species from thermo.dat file and calculate the specific heat, entropy, enthalpy of al the 53 species.

2.To write a program that calculate molecular weights of all the 53 species.

3.To plot the specific heat, entropy, enthaply vs temperature respectively for each species.

MAIN PROGRAM:

close all
clear all
clc

% Inputs
r = 8.314;

% Reading the Thermo.dat file
f1 = fopen('THERMO.dat','r');


%reading header 
file_header = fgetl(f1);


% Getting the global maximum, minimum, mid temperature
temp_range=fgetl(f1);
Global_temp=strsplit(temp_range,' ');
Global_min_temp=str2num(Global_temp{2});
Global_mid_temp=str2num(Global_temp{3});
Global_max_temp=str2num(Global_temp{4});


% Defining working temperature range
t=linspace(Global_min_temp,Global_max_temp,500);


%skipping the three comments
for i=1:3
    skip=fgetl(f1);
end


%running loop for each species there are 53 species in the file
for i=1:53 
    
    %reading the first line that tells which species it is
    %working temp range & no of molecules
    line0=fgetl(f1);
    
    %reading the name of the species
    name=line0(1:6);
    
    %reading the three lines which contains 7 high temp coeff & 7 low temp
    %coeff
    line1=fgetl(f1);
    line2=fgetl(f1);
    line3=fgetl(f1);
    
    %finding the position of E so that we can locate the coefficients
    a=strfind(line1,'E');
    b=strfind(line2,'E');
    c=strfind(line3,'E');
    
    % Extracting Low temp coefficients
    l1=str2num(line1(1:a(1)+3));
    l2=str2num(line1(a(1)+3:a(2)+3));
    l3=str2num(line1(a(2)+4:a(3)+3));
    l4=str2num(line1(a(3)+4:a(4)+3));
    l5=str2num(line1(a(4)+4:a(5)+3));
    l6=str2num(line2(1:b(1)+3));
    l7=str2num(line2(b(1)+4:b(2)+3));
    
    % Extracting High temp coefficients
    h1=str2num(line2(b(2)+4:b(3)+3));
    h2=str2num(line2(b(3)+4:b(4)+3));
    h3=str2num(line2(b(4)+4:b(5)+3));
    h4=str2num(line3(1:c(1)+3));
    h5=str2num(line3(c(1)+4:c(2)+3));
    h6=str2num(line3(c(2)+4:c(3)+3));
    h7=str2num(line3(c(3)+4:c(4)+3));
    
    % Calculating specific heat
    cp=specific_heat(l1,l2,l3,l4,l5,h1,h2,h3,h4,h5,t,r,Global_mid_temp);
    
    %calculating enthalpy
    h=enthalpy(l1,l2,l3,l4,l5,l6,h1,h2,h3,h4,h5,h6,t,r,Global_mid_temp);
    
    %calculating entropy
    s=entropy(l1,l2,l3,l4,l5,l7,h1,h2,h3,h4,h5,h7,t,r,Global_mid_temp);
    
    %plotting
    plotting(cp,h,s,t,name)
    
    %calculating molecular weight of each species
    molecular_weight(i) = mol_wt(name);
end

1.Universal Gas Constant (r) has been given as input that is 8.314 J/mol.k. 2.Opening the thermo.dat file using fopen command, we are reading here that's why we have used 'r' as an permission.

3.Reading first line that is file header using fgetl.

4.Next line used for getting temp_range values and then splitiing into cell array after that we are using these temperature values as Global mid, min max temperatures. Using str2num we are actually converting string value in number if string is there.

5. Defining the working temperature range values using linspace command. 6.Running the for loops which will be calculating Specific Heat, Entropy, Enthalpy, Molecular weight and plotting the graphs for all the 53 species.

7.Reading the first line of every species using fgetl command, then reading the species name here the largest possible species name is upto 6 letter so another fgetl command is used.

8.Reading first, second, third line of every species then after finding the position of 'E' in every line and predicing the lenth of every coefficients.

9.Seperating the coefficients into low and high temperature coefficients. A calling function is made here for Specific Heat, Entropy, Enthalpy, Plotting, Molecular weight.

FUNCTION FOR SPECIFIC HEAT:

function cp=specific_heat(l1,l2,l3,l4,l5,h1,h2,h3,h4,h5,t,r,Global_mid_temp)
if t<=Global_mid_temp
    cp=r*(l1+l2.*(t)+l3.*(t.^2)+l4.*(t.^3)+l5.*(t.^4));
else
    cp=r*(h1+h2.*(t)+h3.*(t.^2)+h4.*(t.^3)+h5.*(t.^4));
end

For this two Specific Heat value for 53 species had been calculated one for high temperature coefficients and one for low temperature coefficients.

FUNCTION FOR ENTROPY:

function s=entropy(l1,l2,l3,l4,l5,l7,h1,h2,h3,h4,h5,h7,t,r,Global_mid_temp)
if t<=Global_mid_temp
    s=r*(l1.*log(t)+l2.*t+l3.*(t.^2)./2+l4.*(t.^3)./3+l5.*(t.^4)./4+l7);
else
    s=r*(h1.*log(t)+h2.*t+h3.*(t.^2)./2+h4.*(t.^3)./3+h5.*(t.^4)./4+h7);
end

For this two Entropy value for 53 species had been calculated one for high temperature coefficients and one for low temperature coefficients.

FUNCTION FOR ENTHALPY:

function h=enthalpy(l1,l2,l3,l4,l5,l6,h1,h2,h3,h4,h5,h6,t,r,Global_mid_temp)
if t<=Global_mid_temp
    h=r*(l1.*t+l2.*(t.^2)./2+l3.*(t.^3)./3+l4.*(t.^4)./4+l5.*(t.^5)./5+l6);
else
    h=r*(h1.*t+h2.*(t.^2)./2+h3.*(t.^3)./3+h4.*(t.^4)./4+h5.*(t.^5)./5+h6);
end

For this two Enthalpy value for 53 species had been calculated one for high temperature coefficients and one for low temperature coefficients.

FUNCTION FOR MOLECULAR WEIGHTS:

function [MW] = mol_wt(species_name)
    Element = {'H' 'C' 'N' 'O' 'R'};
    fprintf('Molecular Weight of ');
    fprintf(species_name);
    fprintf(' : ');
    
    %atomic weight of each individuals
    Element_at_wt = [1 12 14 16 40];
    MW = 0;
    for i = 1:length(species_name)
        for j = 1:length(Element)
            if strcmp(species_name(i), Element(j))
                MW = MW + Element_at_wt(j);
                %Storing the position of matched element
                position = j;
            end
        end
        % This will be used if any string is available in species_name
        n = str2num(species_name(i));
        if n > 1
            MW = MW + Element_at_wt(position)*(n-1);
        end
    end
    fprintf('%d',MW);
    disp(' ');
end

1.All the species is made of elements which are Hydrogen, Carbon, nitrogen, oxygen, Ar. So for these their molecular weights had been taken as input.

2.Starting with for loop one for species name and oher for elements, strcmp command is used to compare the elements with species if both matches then their molecular weight will be calculated alternately with every matching elements with to species.

3.Finding string in species name and converting into number then taking their position for extracting their weight, if 'n' it will be greater than one then it is used in further step to calculate the molecular weight. this 'n' here represents the number of molecules of a element in species.

FUNCTION FOR PLOTTING:

function plotting(cp,h,s,t,species_name)
    cd(pwd)                % Making sure to be in current directory
    dir_name=species_name; % Assigning the folder name
    mkdir(dir_name);       % Creates the folder
    
    % Plotting for Specific Heat and Temperature
    figure(1)
    plot(t,cp,'linewidth',3,'color','r')
    % Using sprintf to print string species name and string in between them
    title(sprintf('Specific Heat vs Temperature for %s:',species_name))
    xlabel('Temperature')
    ylabel('Specific Heat')
    grid on;
    saveas(1,'cp_vs_t.png')
    
    % Plotting for Temperature and Enthalpy
    figure(2)
    plot(t,h,'linewidth',3,'color','b')
    % Using sprintf to print string species name and string in between them
    title(sprintf('Enthalpy vs Temperature for %s:',species_name))
    xlabel('Temperature')
    ylabel('Enthalpy')
    grid on;
    
    %Plotting for Temperature and Entropy    
    figure(3)
    plot(t,s,'linewidth',3,'color','g')
    % Using sprintf to print string species name and string in between them
    title(sprintf('Entropy vs Temperature for %s:',species_name))
    xlabel('Temperature')
    ylabel('Entropy')
    grid on;
    
    % Saving the image in the folder
    cd(dir_name)
    saveas(1,'cp_vs_t.png')
    saveas(2,'h_vs_t.png')
    saveas(3,'s_vs_t.png')
    
    cd ..
end

1.cd(pwd0 used to be in a current directory. 2.Assigning the directory name as species name for every 53 species, and mkdir used for creating for 53 species. As can be seen in the code plot,title,xlabel,ylabel used for every specific heat, entropy, enthalpy vs temperature respectively has been used.

3.Here sprintf used if any string present in the name of species it wood be printed in title. saving all the images as .png file.

NOTE : For everytime we use fgetl command the next line prior to the first line is taken.

PLOT FOR O2:

PLOT FOR N2:

PLOT FOR CO2:

MOLECULAR WEIGHT OF SPECIES: