Project 1 - Parsing NASA thermodynamic data

AIM:-

To write a code in MATLAB/Octave to parse the thermodynamic data file and then calculate the thermodynamic properties of various gas species.

OBJECTIVE:-

1. Write a function that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the data file. You will be reading this file NASA thermodynamic data (Tip: use fopen and fgetl). The formulae are shown below.

Here R is the universal gas constant, T is the temeprature

The FIRST 7 coefficients are HIGH-temperature coefficients and the SECOND 7 coefficients are LOW-temperature coefficients.

Also use the LOCAL TEMPERATURE RANGES which is unique for each individual species, not the GLOBAL TEMPERATURE RANGE.

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

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

4. Save the plots as images with appropriate names and dump them in separate folders for each species with suitable name. All these should be generated automatically.

INTRODUCTION:-

File parsing is the process of analyzing a string of symbols, either in natural language,computer languages or data structures, confroming to the rules of formal grammar . Within computational linguistics the term is used to refer to the formal analysis by a computer of a sentence or other strings of words into its constituents results in a parse tree showing their syntactic relation to each other, which may also contain semantic and other information.

In the following program we are using NASA polynomial format for CHEMKIN-II. Ther chemkin 'Thermo.dat' file format is based very closely on the NASA(1971) file format.

Entropy:

Entropy is a measure of the randomness or disorder of a system.The value of entropy depends on the mass of a system. It is denoted by the letter S and has units of joules per kelvin.

Enthalpy:

Enthalpy is a thermodynamic property of a system. It is the sum of the internal energy added to the product of the pressure and volume of the system. It reflects the capacity to do non-mechanical work and the capacity to release heat.

Specific Heat:

The specific heat is the amount of heat per unit mass required to raise the temperature by one degree Celsius. The relationship between heat and temperature change is usually expressed in the form shown below where c is the specific heat. The relationship does not apply if a phase change is encountered, because the heat added or removed during a phase change does not change the temperature.

Main Program:-

clear all
close all
clc

% Reading the file data 
thermo_data = fopen('THERMO.dat','r');

%universal gas constant 
R = 8.314;

%to read the Header data
fgetl(thermo_data)

%to read the global temperature
global_temp_ran=fgetl(thermo_data)
%To separate the global temperature values
a=strsplit(global_temp_ran,' ');
global_low_temp=str2num(a{1});
global_mid_temp=str2num(a{2});
global_high_temp=str2num(a{3});

temp =linspace(global_low_temp,global_high_temp,400)

%for skipping the next three lines in the data
fgetl(thermo_data);
fgetl(thermo_data);
fgetl(thermo_data);

%To find out species information
for i=1:53
    line1=fgetl(thermo_data)
    A=strsplit(line1,' ');
    species_Name=(A{1});
    
    t2= fgetl(thermo_data);
    t3= fgetl(thermo_data);
    t4= fgetl(thermo_data);
    
    j= strfind(t2,'E');
    k=strfind(t3,'E');
    l=strfind(t4,'E');
    
    %High Temeperature coefficient
    
    h1= str2num(t2(1:j(1)+3));
    h2= str2num(t2(j(1)+4:j(2)+3));
    h3= str2num(t2(j(2)+4:j(3)+3));
    h4= str2num(t2(j(3)+4:j(4)+3));
    h5= str2num(t2(j(4)+4:j(5)+3));
    h6= str2num(t2(1:k(1)+3));
    h7= str2num(t3(k(1)+4:k(2)+3));
    
    % Low temperature coefficient
    n1=str2num(t3(k(2)+4:k(3)+3));
    n2=str2num(t3(k(3)+4:k(4)+3));
    n3=str2num(t3(k(4)+4:k(5)+3));
    n4=str2num(t4(1:l(1)+3));
    n5=str2num(t4(l(1)+4:l(2)+3));
    n6=str2num(t4(l(2)+4:l(3)+3));
    n7=str2num(t4(l(3)+4:l(4)+3));
    
    [c_p]=specific_heat(h1,h2,h3,h4,h5,n1,n2,n3,n4,n5,R,temp,global_mid_temp);
    
    [h_rt]= enthalpy(h1,h2,h3,h4,h5,h6,n1,n2,n3,n4,n5,n6,R,temp,global_mid_temp);
    
    [S]= entropy(h1,h2,h3,h4,h5,h6,h7,n1,n2,n3,n4,n5,n6,n7,R,temp,global_mid_temp);
    
    mkdir(species_Name);
    current_dir=pwd;
    cd(species_Name)
    
    figure(1)
    plot(temp,c_p,'linewidth',2,'color','g');
    title(species_Name)
    xlabel("Temperature")
    ylabel("Specific Heat")
    test =sprintf('frame_cp%d.png',i);
    saveas (gca,test);
    
    figure(2)
    plot(temp,h_rt,'linewidth',2,'color','b');
    title(species_Name)
    xlabel('Temperature')
    ylabel("Enthapy")
    test =sprintf('frame_h%d.png',i);
    saveas (gca,test);
    
    figure(3)
    plot(temp,S,'linewidth',2,'color','r');
    title(species_Name)
    xlabel('Temperature');
    ylabel('Entropy');
    test =sprintf('frame_s%d.png',i);
    saveas (gca,test);
    
    mw = molecular_weight(species_Name);
    
    cd(current_dir);
    
    
end

Function code for the calcuation of specific heat:-

function[c_p]=specific_heat(h1,h2,h3,h4,h5,n1,n2,n3,n4,n5,R,temp,global_mid_temp);

for i =1:400
  if temp(i) <=global_mid_temp
  
  c_p(i)=R*(n1+(n2*temp(i))+(n3*(temp(i)).^2)+(n4*(temp(i)).^3)+(n5*(temp(i)).^4));
else
  c_p(i)=R*(h1+(h2*temp(i))+(h3*(temp(i)).^2)+(h4*(temp(i)).^3)+(h5*(temp(i)).^4));
endif
endfor
end

Function code for the calculation of Entropy:-

function [S] = entropy(h1,h2,h3,h4,h5,h6,h7,l1,l2,l3,l4,l5,l6,l7,R,temp,global_mid_temp)
  for i= 1:400
    if temp(i)<=global_mid_temp
    S(i)=R.*(l1.*log(temp(i)))+(l2.*temp(i))+(l3.*(((temp(i)).^2)./2))+(l4.*(((temp(i)).^3)./3))+(l5.*(((temp(i)).^4))+l7);
  else
    S(i)=R.*(h1.*log(temp(i)))+(h2.*temp(i))+(h3.*(((temp(i)).^2)./2))+(h4.*(((temp(i)).^3)./3))+(h5.*(((temp(i)).^4))+h7);
    end
  endfor
endfunction

Function code for the calculation of Enthalpy:-

function[h_rt]=enthalpy(h1,h2,h3,h4,h5,h6,n1,n2,n3,n4,n5,n6,R,temp,global_mid_temp)

for i=1:400
  if temp(i) <=global_mid_temp
   
   h_rt(i)=R.*temp(i).*(n1+(n2.*temp(i)./2))+(n3.*(temp(i)./3))+(n4.*(temp(i)./4))+(n5.*(temp(i)./5))+(n6./temp(i));
 else
   h_rt(i)=R.*temp(i).*h1+(h2.*(temp(i)./2))+(h3.*(temp(i)./3))+(h4.*(temp(i)./4))+(h5.*(temp(i)./5))+(h6./temp(i));
end
end
end

Function for the calculation of Molecular Weight:-

function[mw]=molecular_weight(species_Name)

elements={'H','C','O','N','A','S'};
atomic_wt=[1 12.01 15.99 14 39.95 32.06];
mw=0;
for i=1:length(species_Name);
  for j=1:length(elements);
    if strcmpi(species_Name(i),elements(j))
      mw= mw+atomic_wt(j)
      q=j;
    end  
    endfor
    m=str2num(species_Name(i));
    if m>1
      mw=mw+(atomic_wt(q)*(m-1));
     
    end
endfor

Errors Encountered:-

• Syntax errors.
• Errors where encountered while separating variables name.
• Erros where while creating the folder for the storage of output.

OUTPUT:-

plots of CH₂CHO

screenshot to show location of folders created for species.