Challenge

Project 1 - Parsing NASA thermodynamic data

MMd.Rubayat Alam Riad28 Jan 2021

AIM:- MATLAB program to fulfill the following objectives :

Writing a code in MATLAB to parse the NASA thermodynamic data file.
Writing a function that extracts the 14 coefficients and calculates the enthalpy, entropy, and specific heats for all the species in the data file.
Calculating the molecular weight of each species and displaying it in the command window.
Plotting the Cp, Enthalpy, and Entropy for the local temperature range (low temperature: high temperature) specific for each species.
Saving the plots as images with appropriate names and dumping them in separate folders for each species with a suitable name.

MAIN CONCEPT: 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.

The thermodynamic properties in terms of the coefficients are given as :

The first 7 coefficients are High – temperature coefficients and the second 7 coefficients are low- temperature coefficients.

Here R is the universal gas constant, T is the temperature.

Specific Heat: It is defined as the amount of heat required to raise the temperature of a unit mass of a substance by 1 degree.

Enthalpy: It is the sum of the internal energy and the product of the pressure and volume of a thermodynamic process. Enthalpy is an energy-like property or state function—it has the dimensions of energy (and is thus measured in units of joules or ergs), and its value is determined entirely by the temperature, pressure, and composition of the system and not by its history.

Entropy: It is 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. The concept of entropy provides deep insight into the direction of spontaneous change for many everyday phenomena.

File parsing here means extracting useful information from a file and utilizing that information for our program.

SOLUTION ALGORITHM

The program begins with the clear all, close all, and clc commands to free up the command window and the workspace.
We will then read the data from the external NASA thermodynamic file using fopen(to open the file) and fgetl(to read the content in the file line by line ) commands.
The useful header information will be extracted from the file such as global temperature limits.
Then we will proceed to extract the individual species data of all the species present in the file using fgetl command to read the file line by line and skipping the line not containing any useful data.
The data that is to be extracted for all the 53 species mentioned in the file are 14 coefficients(7 for higher temperature and & for lower temperature), species name, and the local temperature limits.
This process will be carried out with the help of a for aloop.
Calculation of specific heat, enthalpy, entropy, and molecular weight will be done in separate function files and will be called in the main program.
Then we will create plots for specific heat v/s temperature, enthalpy v/s temperature, and entropy v/s temperature for all the 53 species and these 3 plots will be saved in separate folders for each species.
After this, we will display the molecular weight of all the species in the command window.
The code for the entire program will be written in MATLAB.

MATLAB CODE

THE ROOT PROGRAM:

%program for file parsing of NASA thermodynamic file
close all
clear all
clc

%opening the external file
f1= fopen('THERMO.dat','r')
line_1= fgetl(f1); %skipping this line
line_2= fgetl(f1);

%obtaining global temperature values 
%storing them as numbers in seperate values 
global_values= strsplit(line_2,' ');
g_low_temp= str2num(global_values{2})
g_mid_temp= str2num(global_values{3})
g_high_temp= str2num(global_values{4})

line_3= fgetl(f1);%skipping the line
line_4= fgetl(f1);%skipping the line
line_5= fgetl(f1);

%defining input universal gas constant in KJ/K-mol-K 
R= 8.314;

for i= 1:53
    T= fgetl(f1);
    
    %extracting the species name
    species_name= T(1:9);
    species_name= deblank(species_name); %removing blank spaces
    
    %extracting local temperatures
    species_temp= T(49:75);
    local_values= strsplit(species_temp,' ');
    l_low_temp= str2num(local_values{1});
    l_mid_temp= str2num(local_values{3});
    l_high_temp= str2num(local_values{2});
    %storing the local temperatures in a variable
    temp= [l_low_temp l_mid_temp l_high_temp];
    
    temp_range=linspace(l_low_temp,l_high_temp,500);
    
    coeff= fgetl(f1);
    coeff1= findstr(coeff,'E');
    
    %extracting the coefficients
    a1= coeff(1:coeff1(1)+3);
    a1= str2num(a1);
    a2= coeff(coeff1(1)+4:coeff1(2)+3);
    a2= str2num(a2);
    a3= coeff(coeff1(2)+4:coeff1(3)+3);
    a3= str2num(a3);
    a4= coeff(coeff1(3)+4:coeff1(4)+3);
    a4= str2num(a4);
    a5= coeff(coeff1(4)+4:coeff1(5)+3);
    a5= str2num(a5);
    
    new_line= fgetl(f1);
    coeff2= findstr(new_line,'E');
    a6= new_line(1:coeff2(1)+3);
    a6= str2num(a6);
    a7= new_line(coeff2(1)+4:coeff2(2)+3);
    a7= str2num(a7);
    a8= new_line(coeff2(2)+4:coeff2(3)+3);
    a8= str2num(a8);
    a9= new_line(coeff2(3)+4:coeff2(4)+3);
    a9= str2num(a9);
    a10= new_line(coeff2(4)+4:coeff2(5)+3);
    a10= str2num(a10);
    
   new_line1= fgetl(f1);
   coeff_3= findstr(new_line1,'E');
   a11= new_line1(1:coeff_3(1)+3);
   a11= str2num(a11);
   a12= new_line1(coeff_3(1)+4:coeff_3(2)+3);
   a12= str2num(a12);
   a13= new_line1(coeff_3(2)+4:coeff_3(3)+3);
   a13= str2num(a13);
   a14= new_line1(coeff_3(3)+4:coeff_3(4)+3);
   a14= str2num(a14);
   
  %specific heat calculation
  Cp= specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,temp_range,R);
   
  %enthalpy calculation
  H= enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a2,a13,temp,temp_range,R);
  
  %entropy calculation
  S= entropy_1(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,temp_range,R);
  
  %calculation of molecular weight
  m_weight= weight(species_name);
  
  %saving molecular weight
  species_w= fopen('species_mass.txt','w');
  fprintf(species_w,'molecular mass of %s is %f',species_name,m_weight);
  fclose(species_w);
  
%plotting and saving
cur_dir=pwd;

%making a new folder for each species
mkdir(species_name);
cd(species_name);


figure(1)
plot(temp_range,Cp,'linewidth',2,'color',[0.2345 0.4567 0.5678])
xlabel('Temperature(K)');
ylabel('Specific heat(KJ/kg-K)');
title(species_name);

figure(2), 
plot(temp_range,H,'LineWidth',2,'color',[0.5678 0.3423 0.7874])
xlabel('Temperature(K)');
ylabel('Enthalpy(KJ/Kg)');
title(species_name);

figure(3)
plot(temp_range,S,'LineWidth',2,'color',[0.4532 0.1234 0.3324])
xlabel('Temperature(K)');
ylabel('Entropy(KJ/K)');
title(species_name);


%saving the plots in spearate folder for each species
saveas(1,'specific heat.jpg');
saveas(2,'enthalpy.jpg');
saveas(3,'entropy.jpg');

cd(cur_dir);

%printing the molcular weight of the species
sprintf('molecular weight of %s = %f',species_name,m_weight)


end

%closing the file
fclose(f1)

Enthalpy calculation cod:

function[S]= entropy_1(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,temp_range,R)
    
    %a1 a2 a3 a4 a5 are high temperature coefficients
    %a8 a9 a10 a11 a12 are low temperature coefficients
    
    for i= 1:length(temp_range)
        if temp_range(i)<=temp(2)
            S(i)= R.*((a8*log(temp_range(i))) + (a9*temp_range(i)) + (a10*((temp_range(i).^2)/2)) + (a11*((temp_range(i).^3)/3)) + (a12*((temp_range(i).^4)/4)) + a14);
            
        else
            S(i)= R.*((a1*log(temp_range(i))) + (a2*temp_range(i)) + (a3*((temp_range(i).^2)/2)) + (a4*((temp_range(i).^3)/3)) + (a5*((temp_range(i).^4)/4)) + a7);
            
        end
end

Entropy calculation code:

function[S]= entropy_1(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,temp_range,R)
    
    %a1 a2 a3 a4 a5 are high temperature coefficients
    %a8 a9 a10 a11 a12 are low temperature coefficients
    
    for i= 1:length(temp_range)
        if temp_range(i)<=temp(2)
            S(i)= R.*((a8*log(temp_range(i))) + (a9*temp_range(i)) + (a10*((temp_range(i).^2)/2)) + (a11*((temp_range(i).^3)/3)) + (a12*((temp_range(i).^4)/4)) + a14);
            
        else
            S(i)= R.*((a1*log(temp_range(i))) + (a2*temp_range(i)) + (a3*((temp_range(i).^2)/2)) + (a4*((temp_range(i).^3)/3)) + (a5*((temp_range(i).^4)/4)) + a7);
            
        end
end

Molecular weight calculation code:

function[m_weight]= weight(species_name)
    
    element= {'H','O','N','Ar','C','S'};
    atomic_mass= [1.008, 15.999, 14.007, 39.948, 12.011, 32.065];
    m_weight=0;
    for i= 1:length(species_name)
        for j= 1:length(element)
            if strcmp(species_name(i),element(j))
                m_weight= m_weight + atomic_mass(j);
                a= j;
            end
                
            end
            
        b= str2num(species_name(i));
        if b>1
            m_weight= m_weight + atomic_mass(a)*(b-1);
        end
    end
end

OUTPUT

After running the code the following output will be obtained:

f1 =

     3


g_low_temp =

   300


g_mid_temp =

        1000


g_high_temp =

        5000


ans =

    'molecular weight of O = 15.999000'


ans =

    'molecular weight of O2 = 31.998000'


ans =

    'molecular weight of H = 1.008000'


ans =

    'molecular weight of H2 = 2.016000'


ans =

    'molecular weight of OH = 17.007000'


ans =

    'molecular weight of H2O = 18.015000'


ans =

    'molecular weight of HO2 = 33.006000'


ans =

    'molecular weight of H2O2 = 34.014000'


ans =

    'molecular weight of C = 12.011000'


ans =

    'molecular weight of CH = 13.019000'


ans =

    'molecular weight of CH2 = 14.027000'


ans =

    'molecular weight of CH2(S) = 46.092000'


ans =

    'molecular weight of CH3 = 15.035000'


ans =

    'molecular weight of CH4 = 16.043000'


ans =

    'molecular weight of CO = 28.010000'


ans =

    'molecular weight of CO2 = 44.009000'


ans =

    'molecular weight of HCO = 29.018000'


ans =

    'molecular weight of CH2O = 30.026000'


ans =

    'molecular weight of CH2OH = 31.034000'


ans =

    'molecular weight of CH3O = 31.034000'


ans =

    'molecular weight of CH3OH = 32.042000'


ans =

    'molecular weight of C2H = 25.030000'


ans =

    'molecular weight of C2H2 = 26.038000'


ans =

    'molecular weight of C2H3 = 27.046000'


ans =

    'molecular weight of C2H4 = 28.054000'


ans =

    'molecular weight of C2H5 = 29.062000'


ans =

    'molecular weight of C2H6 = 30.070000'


ans =

    'molecular weight of CH2CO = 42.037000'


ans =

    'molecular weight of HCCO = 41.029000'


ans =

    'molecular weight of HCCOH = 42.037000'


ans =

    'molecular weight of H2CN = 28.034000'


ans =

    'molecular weight of HCN = 27.026000'


ans =

    'molecular weight of HNO = 31.014000'


ans =

    'molecular weight of N = 14.007000'


ans =

    'molecular weight of NNH = 29.022000'


ans =

    'molecular weight of N2O = 44.013000'


ans =

    'molecular weight of NH = 15.015000'


ans =

    'molecular weight of NH2 = 16.023000'


ans =

    'molecular weight of NH3 = 17.031000'


ans =

    'molecular weight of NO = 30.006000'


ans =

    'molecular weight of NO2 = 46.005000'


ans =

    'molecular weight of HCNO = 43.025000'


ans =

    'molecular weight of HOCN = 43.025000'


ans =

    'molecular weight of HNCO = 43.025000'


ans =

    'molecular weight of NCO = 42.017000'


ans =

    'molecular weight of CN = 26.018000'


ans =

    'molecular weight of HCNN = 41.033000'


ans =

    'molecular weight of N2 = 28.014000'


ans =

    'molecular weight of AR = 0.000000'


ans =

    'molecular weight of C3H8 = 44.097000'


ans =

    'molecular weight of C3H7 = 43.089000'


ans =

    'molecular weight of CH3CHO = 44.053000'


ans =

    'molecular weight of CH2CHO = 43.045000'


ans =

     0

>>

Ploting for O2 , N2 & CO2 are shown bellows:

For O2