Project 1 - Parsing NASA thermodynamic data

Parsing NASA thermodynamic data

Objective :

1. To write a function that extracts the 14 coefficients and calculates the enthalpy, entropy and specific heats for all the species in the NASA thermodynamic data file "THERNO.dat".

The NASA polynomials have the form as shown below:

Where R = universal gas constant and T = Temperature

2. To calculate the molecular weight of each species and display it in the command window.

3. To plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species.

4. To save the plots as images with appropriate names and dump them in separate folders for each species with suitable name automatically.

File parsing :

File parsing refers to the conversion of data from the external file to a soft data format that can be used in a computer program.

NASA thermodynamic file :

The given THERMO.dat file contains 14 coefficients for 53 different elements. These coefficients are used to calculate the specific heat, enthalpy and entropy for each element. Each element has its own range of local temperature. The parent elements contained in the file are Argon, Oxygen, Hydrogen, Nitrogen and Carbon.

For analyzing a thermodynamic processes, we must need enthalpy, entropy and specific heats of various species. For example, in a combustion process, we required to observe that how the combustion is going on and how efficient it is. It observe this, we required to study the thermodynamic properties of various species that are engaged in combustion.

NASA came up with polynomials that can be used to evaluate thermodynamic properties. They have documented the coefficients that are required to evaluate these polynomials for more than 1000 species.

The NASA polynomials have the form :

In this, there are total 14 coefficients are present, first 7 are high-temperature coefficients (>1000) and second 7 are low-temperature coefficients (<1000).

CODE :

1. Specific heat :

The amount of energy that must be added in the form of heat to the one unit of mass of the substance in order increase the temperature is called specific heat.

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

Specific heat function code :

% ========================= Specific Heat function file =========================== %

function Cp = specificheat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp)

% comparing with global temperature %
% If less than global temp, min coefficents are taken %
% otherwise max coefficients are taken %

if Temp <= global_medium_temp
    Cp = R*(a8+(a9*Temp)+(a10*(Temp).^2)+(a11*(Temp).^3)+(a12*(Temp).^4));
else
    Cp = R*(a1+(a2*Temp)+(a3*(Temp).^2)+(a4*(Temp).^3)+(a5*(Temp).^4));
end

2. Enthalpy :

It is a property of a thermodynamic system, is equal to the sum of internal energy and the product of its pressure and volume. The heat absorbed or released is equal to the change in internal energy.

H = (a1 + a2*T/2 + a3*T^2/3 + a4*T^3/4 + a5*T^4/5 +a6/T)*R*T

Enthalpy function code :

% ======================================= Enthalpy function file ======================================== %
function H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,Temp,R,global_medium_temp)

% comparing with global temperature %
% If less than global temp, min coefficents are taken %
% otherwise max coefficients are taken %

if Temp <= global_medium_temp
    H = R*Temp.*(a8+((a9*Temp)/2)+(a10*((Temp).^2)/3)+(a11*((Temp).^3)/4)+(a12*((Temp).^4)/5)+(a13./Temp));
else
    H = R*Temp.*(a1+((a2.*Temp)/2)+(a3*((Temp).^2)/3)+(a4*((Temp).^3)/4)+(a5*((Temp).^4)/5)+(a6./Temp));
end

3. Entropy :

It is defined as the amount the energy that would be release per unit temperature which is unavailable for the useful work.

S = R*(a1 ln*T + a2*T +a3*T^2/2 + a4*T^3/3 + a5*T^4/4 + a7)

Entropy function code :

% ================================== Entropy function file ======================================= %
function S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,Temp,R,global_medium_temp)

% comparing with global temperature %
% If less than global temp, min coefficents are taken %
% otherwise max coefficients are taken %

if Temp <= global_medium_temp
    S = R*(a8.*log(Temp)+(a9*Temp)+(a10*((Temp).^2)/2)+(a11*((Temp).^3)/3)+(a12*((Temp).^4)/4)+a14);
else
    S = R*(a1.*log(Temp)+(a2.*Temp)+(a3*((Temp).^2)/2)+(a4*((Temp).^3)/3)+(a5*((Temp).^4)/4)+a7);
end

4. Molecular weight :

The mass of molecule of the substance is called as molecular mass.

Molecular weight function code :

% Molecular weight function file %
function molecularweight = molecular_weight(Species_SP);
% standard atomic names - hydrogen,carbon,oxygen,nitrogen,argon
atomic_name = ['H','C','O','N','A'];

% standard atomic numbers
atomic_weight = [1.008 12.011 15.999 14.007 39.948];

% starting range 
molecularweight = 0;

% 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 of %s = %d',Species_SP,molecularweight);
 disp(' ');

end

MAIN CODE :

clear all
close all
clc

% reading given thermo.dat file %
f1 = fopen('THERMO.dat','r');
gt = fgetl(f1);

% reading header information %
global_temp = fgetl(f1);

% splitting the temperature data %
t = strsplit(global_temp);

% converting strings to numbers %
global_low_temp = str2num(t{2});
global_medium_temp = str2num(t{3});
global_high_temp = str2num(t{4});

% temperature range is created b/w low and high temperatures %
Temp = linspace(global_low_temp,global_high_temp,1000);

% universal gas constant
R = 8.314;

% skipping first three lines in the given file which contain comments %
for i = 1:3
    gt = fgetl(f1);
end

% reading the 53 species %
for i = 1:53
    gt = fgetl(f1);
    A = strsplit(gt,' ');
    Species_SP = (A{1});
    fprintf('%s ', '>>> Element name :','"',Species_SP,'"')
    
    line_1 = fgetl(f1);
    a = strfind(line_1,'E');
    
    % extracting first 7 high temperature coefficents %
    a1 = line_1(1:a(1)+3);
    a1 = str2num(a1);
    a2 = line_1(a(1)+4:a(2)+3);
    a2 = str2num(a2);
    a3 = line_1(a(2)+4:a(3)+3);
    a3 = str2num(a3);
    a4 = line_1(a(3)+4:a(4)+3);
    a4 = str2num(a4);
    a5 = line_1(a(4)+4:a(5)+3);
    a5 = str2num(a5);
    
    line_2 = fgetl(f1);
    b = strfind(line_1,'E');
    a6 = line_2(1:b(1)+3);
    a6 = str2num(a6);
    a7 = line_2(a(1)+4:b(2)+3);
    a7 = str2num(a7);
    
    % extracting second 7 low temperature coefficents %
    a8 = line_2(b(2)+4:b(3)+3);
    a8 = str2num(a8);
    a9 = line_2(b(3)+4:b(4)+3);
    a9 = str2num(a9);
    a10 = line_2(b(4)+4:b(5)+3);
    a10 = str2num(a10);
    
    line_3 = fgetl(f1);
    c = strfind(line_1,'E');
    
    a11 = line_3(1:c(1)+3);
    a11 = str2num(a11);
    a12 = line_3(c(1)+4:c(2)+3);
    a12 = str2num(a12);
    a13 = line_3(c(2)+4:c(3)+3);
    a13 = str2num(a13);
    a14 = line_3(c(3)+4:c(4)+3);
    a14 = str2num(a14);
    
    % Specific heat %
    Cp = specificheat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp);
    
    % Enthalpy %
    H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,Temp,R,global_medium_temp);

    % Entropy %
    S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,Temp,R,global_medium_temp);

    % Molecular weight
    molecularweight = molecular_weight(Species_SP);
        
    % read & write molecular weights %
    molecular_list = fopen('molecular weight.txt','a+');
    fprintf(molecular_list,'molecular weight of %s is %d',Species_SP,molecularweight);
    fclose(molecular_list);
    
    % plotting (Cp,H,S,Temp,Species_SP) and saving as folder %
    cur_dir = pwd;
    mkdir(Species_SP);
    cd(Species_SP);
    
    % plot 1 - Cp vs Temp %
    figure(1)
    plot(Temp,Cp,'linewidth',2,'color','r');
    xlabel('Temperature');
    ylabel('Specific Heat');
    title(sprintf('Variation of Specific Heat with Temperature Range - %s',Species_SP));
    grid on
    saveas(figure(1),'specific heat.jpg');
    
    % plot 2 - E vs Temp %
    figure(2)
    plot(Temp,H,'linewidth',2,'color','b');
    xlabel('Temperature')
    ylabel('Enthalpy')
    title(sprintf('Variation of Enthalpy with Temperature Range - %s',Species_SP));
    grid on
    saveas(figure(2),'enthalpy.jpg')
    
    % plot 3 - S vs Temp %
    figure(3)
    plot(Temp,S,'linewidth',2,'color','g');
    xlabel('Temperature')
    ylabel('Entropy')
    title(sprintf('Variation of Entropy with Temperature Range - %s',Species_SP));
    grid on
    saveas(figure(3),'entropy.jpg')
    
    cd(cur_dir);
    
    % molecuar weight %
    M = molecular_weight(Species_SP);
    f2 = fopen('molecular mass','w');
    fprintf('%s of %s is %f ',' -   Molecular Weight',Species_SP,M)
    fclose(f2)
end

CODE explanation :

1. The code was started with the basic commands clear all, close all and clc.
2. Open the given THERMO.dat file using the 'fopen' command and use 'r' to read the file.
3. 'fgetl' command is used to read the line in the text file.
4. To get the global temperature ranges from the given text file, use the 'strsplit' command to split the temperature into a particular string.
5. Convert the obtained string cell into the numbers using 'str2num' command.
6. A temperature range cell is created between the global low-temperature and global high-temperature using 'linspace' command.
7. Universal gas constant (R) value input given as 8.314.
8. 'For loop' is used to skip the first three lines in the given text file which containing the comments.
9. Again 'for loop' is used to read all 53 species contained in the text file.
10. To find the strings, 'strfind' command is used.
11. All the 14 coefficients are found using the 'strfind' and they are converted into the numbers using 'str2num'.
12. Saved function files of Specific heat, Enthalpy, Entropy and Molecular weight are called to find them.
13. Creating directory using 'mkdir' command.
14. Plots are created between temperature and specific heat, enthalpy and entropy, and they are saved in the folder with their species name in jpg format.
15. Molecular weight is displayed on the command window along with species name using 'fprintf' command.

OUTPUTS :

Molecular weights :

The below is the copy of output obtained from the command window

>>> Element name : " O "  -   Molecular Weight of O is 15.999000 
ans =

     0

>>> Element name : " O2 "  -   Molecular Weight of O2 is 31.998000 
ans =

     0

>>> Element name : " H "  -   Molecular Weight of H is 1.008000 
ans =

     0

>>> Element name : " H2 "  -   Molecular Weight of H2 is 2.016000 
ans =

     0

>>> Element name : " OH "  -   Molecular Weight of OH is 17.007000 
ans =

     0

>>> Element name : " H2O "  -   Molecular Weight of H2O is 18.015000 
ans =

     0

>>> Element name : " HO2 "  -   Molecular Weight of HO2 is 33.006000 
ans =

     0

>>> Element name : " H2O2 "  -   Molecular Weight of H2O2 is 34.014000 
ans =

     0

>>> Element name : " C "  -   Molecular Weight of C is 12.011000 
ans =

     0

>>> Element name : " CH "  -   Molecular Weight of CH is 13.019000 
ans =

     0

>>> Element name : " CH2 "  -   Molecular Weight of CH2 is 14.027000 
ans =

     0

>>> Element name : " CH2(S) "  -   Molecular Weight of CH2(S) is 14.027000 
ans =

     0

>>> Element name : " CH3 "  -   Molecular Weight of CH3 is 15.035000 
ans =

     0

>>> Element name : " CH4 "  -   Molecular Weight of CH4 is 16.043000 
ans =

     0

>>> Element name : " CO "  -   Molecular Weight of CO is 28.010000 
ans =

     0

>>> Element name : " CO2 "  -   Molecular Weight of CO2 is 44.009000 
ans =

     0

>>> Element name : " HCO "  -   Molecular Weight of HCO is 29.018000 
ans =

     0

>>> Element name : " CH2O "  -   Molecular Weight of CH2O is 30.026000 
ans =

     0

>>> Element name : " CH2OH "  -   Molecular Weight of CH2OH is 31.034000 
ans =

     0

>>> Element name : " CH3O "  -   Molecular Weight of CH3O is 31.034000 
ans =

     0

>>> Element name : " CH3OH "  -   Molecular Weight of CH3OH is 32.042000 
ans =

     0

>>> Element name : " C2H "  -   Molecular Weight of C2H is 25.030000 
ans =

     0

>>> Element name : " C2H2 "  -   Molecular Weight of C2H2 is 26.038000 
ans =

     0

>>> Element name : " C2H3 "  -   Molecular Weight of C2H3 is 27.046000 
ans =

     0

>>> Element name : " C2H4 "  -   Molecular Weight of C2H4 is 28.054000 
ans =

     0

>>> Element name : " C2H5 "  -   Molecular Weight of C2H5 is 29.062000 
ans =

     0

>>> Element name : " C2H6 "  -   Molecular Weight of C2H6 is 30.070000 
ans =

     0

>>> Element name : " CH2CO "  -   Molecular Weight of CH2CO is 42.037000 
ans =

     0

>>> Element name : " HCCO "  -   Molecular Weight of HCCO is 41.029000 
ans =

     0

>>> Element name : " HCCOH "  -   Molecular Weight of HCCOH is 42.037000 
ans =

     0

>>> Element name : " H2CN "  -   Molecular Weight of H2CN is 28.034000 
ans =

     0

>>> Element name : " HCN "  -   Molecular Weight of HCN is 27.026000 
ans =

     0

>>> Element name : " HNO "  -   Molecular Weight of HNO is 31.014000 
ans =

     0

>>> Element name : " N "  -   Molecular Weight of N is 14.007000 
ans =

     0

>>> Element name : " NNH "  -   Molecular Weight of NNH is 29.022000 
ans =

     0

>>> Element name : " N2O "  -   Molecular Weight of N2O is 44.013000 
ans =

     0

>>> Element name : " NH "  -   Molecular Weight of NH is 15.015000 
ans =

     0

>>> Element name : " NH2 "  -   Molecular Weight of NH2 is 16.023000 
ans =

     0

>>> Element name : " NH3 "  -   Molecular Weight of NH3 is 17.031000 
ans =

     0

>>> Element name : " NO "  -   Molecular Weight of NO is 30.006000 
ans =

     0

>>> Element name : " NO2 "  -   Molecular Weight of NO2 is 46.005000 
ans =

     0

>>> Element name : " HCNO "  -   Molecular Weight of HCNO is 43.025000 
ans =

     0

>>> Element name : " HOCN "  -   Molecular Weight of HOCN is 43.025000 
ans =

     0

>>> Element name : " HNCO "  -   Molecular Weight of HNCO is 43.025000 
ans =

     0

>>> Element name : " NCO "  -   Molecular Weight of NCO is 42.017000 
ans =

     0

>>> Element name : " CN "  -   Molecular Weight of CN is 26.018000 
ans =

     0

>>> Element name : " HCNN "  -   Molecular Weight of HCNN is 41.033000 
ans =

     0

>>> Element name : " N2 "  -   Molecular Weight of N2 is 28.014000 
ans =

     0

>>> Element name : " AR "  -   Molecular Weight of AR is 39.948000 
ans =

     0

>>> Element name : " C3H8 "  -   Molecular Weight of C3H8 is 44.097000 
ans =

     0

>>> Element name : " C3H7 "  -   Molecular Weight of C3H7 is 43.089000 
ans =

     0

>>> Element name : " CH3CHO "  -   Molecular Weight of CH3CHO is 44.053000 
ans =

     0

>>> Element name : " CH2CHO "  -   Molecular Weight of CH2CHO is 43.045000 
ans =

     0

>> 

Screenshot of window showing folders :

RESULTS :

1. Plots for O2

2. Plots for N2

3. Plots for CO2

Google Drive link for the output folders :

https://drive.google.com/drive/folders/1DYSWLM0Ls_mv7Rv3s7GEgj7sWp8GVDwG?usp=sharing

CONCLUSION :

From the given NASA thermodynamic data, all the coefficients are found and molecular weights of each element is calculated using MATLAB. Also, plots are created for temperature and specific heat, enthalpy and entropy. All this data is saved in individual folders.

ERRORS :

Error 1 :

The molecular weight is not specified above, so this error was created.

Error 2 :

When I run the code, the output for molecular weight is displayed twice. I did not get that why it was happened. This error was simply created by executing and printing the same code in both function code and main code.

THANK YOU