Project 1 - Parsing NASA thermodynamic data

NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S. They have also documented the co-efficients that are required to evaluate these polynomials for 1000+ species.

To 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.

The formulae are shown below.

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

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

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

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

To 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.

The main program that you will submit will just take the "THERMO.dat" file as an input and generate Cp, H and S plots.

Main program for parsing NASA thermodynamic data

clear all
close all
clc

%Opening of NASA's thermodynamic data

f1 = fopen('THERMO.dat','r');
global_temp = '300.000  1000.000  5000.000';

%Obtainingof globalmax, min and medium temperature
first_line = fgetl(f1);
temp_range = strsplit(global_temp,' ');
global_min_temp = str2num(temp_range{1})
global_med_temp = str2num(temp_range{2})
global_max_temp = str2num(temp_range{3})


% Range of temperature
T = linspace(global_min_temp,global_max_temp,400);
R = 8.314;

% skippingthe next three comments in the file
for i=3:6
    comments = fgetl(f1);
end

% getting the coeffs of species

for i=1:53
 %Extracting line 1 and the name of the species
 line1 = fgetl(f1);
 name = strtok(line1);

 %Extracting higher and lower temperature coeffs
 line2 = fgetl(f1);
 line3 = fgetl(f1);
 line4 = fgetl(f1);

 A=strfind(line2,'E');
 B=strfind(line3,'E'); 
 C=strfind(line4,'E');

 %Extracting higher temperature coeffs
 a1 = str2num(line2(1:A(1)+3));
 a2 = str2num(line2(A(1)+4:A(2)+3));
 a3 = str2num(line2(A(2)+4:A(3)+3));
 a4 = str2num(line2(A(3)+4:A(4)+3));
 a5 = str2num(line2(A(4)+4:A(5)+3));
 a6 = str2num(line3(1:B(1)+3));
 a7 = str2num(line3(B(1)+4:B(2)+3));

 %Extracting higher temperature coeffs
 b1 = str2num(line3(B(2)+4:B(3)+3));
 b2 = str2num(line3(B(3)+4:B(4)+3));
 b3 = str2num(line3(B(4)+4:B(5)+3));
 b4 = str2num(line4(1:C(1)+3));
 b5 = str2num(line4(C(1)+4:C(2)+3));
 b6 = str2num(line4(C(2)+4:C(3)+3));
 b7 = str2num(line4(C(3)+4:C(4)+3));

 %Specific heat calculation 
 cp = specificheat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,global_med_temp);
 
 %Enthalpy calculation
 h = enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T,R,global_med_temp);

 % Entropy calculation
 s = entropy_calc(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T,R,global_med_temp);

 % Calculating molecular weight of species
 Molecular_wt(i) = molecular_weight(name);

 % Calculating individual folder for each species and plotting
 plotting(cp,h,s,T,name)
end

For molecular weight display

function MW = molecular_weight(name)
atoms=['H','C','O','N','R'];
 fprintf('Molecular Weight of  ');
 fprintf(name);
 fprintf(' : ');
 atomic_weight = [1 12 16 14 40];
 MW = 0;
 for i = 1:length(name)
     for j = 1:length(atoms) 
         if strcmp(name(i),atoms(j))
             MW =MW+atomic_weight(j);
             position = j;
         end
     end
      n = str2num(name(1));
     if n>1
         MW = MW + atomic_weight(position)*(n-1);
     end
 end
 fprintf('%d',MW);
 disp(' ');
end

For plotting the graphs of enthalpy,entropy & specific heat for various molecules of atom

% Creating plot  function
function plotting(cp,h,s,T,name)
cd(pwd)
dir_name = name;
mkdir(dir_name);

figure(1)
plot(T,cp,'linewidth',2,'color','b')
title(sprintf('Specific heat Vs Temperature for %s',name))
xlabel('TEMPERATURE')
ylabel('SPECIFIC HEAT')

figure(2)
plot(T,h ,'linewidth',2,'color','r')
title(sprintf('Enthaply Vs Temperature for %s',name))
xlabel('TEMPERATURE')
ylabel('Enthaply')

figure(3)
plot(T,s,'linewidth',2,'color','g')
title(sprintf('Entropy Vs Temperature for %s',name))
xlabel('TEMPERATURE')
ylabel('Entropy')

cd(dir_name)

  saveas(1,'cp_vs_t.png')
  saveas(2,'h_vs_t.png')
  saveas(3,'s_vs_t.png')
   
   cd .. 
end

For entropy calculation

% Calculation of entropy
function s = entropy_calc(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T,R,global_med_temp)
 if T<= global_med_temp
     s = R*(a1.*log(T)+a2.*T+a3.*(T.^2)./2+a4.*(T.^3)./3+a5.*(T.^4)./4+a7);
 else
     s = R*(b1.*log(T)+b2.*T+b3.*(T.^2)./2+b4.*(T.^3)./3+b5.*(T.^4)./4+b7);
 end

For enthalpy calculation

% Calculation of enthaply
function  h= enthapy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T,R,global_med_temp)
 if T<= global_med_temp
     h = R.*(a1.*T+(a2.*T.^2)./2+(a3.*T.^3)./3+(a4.*T.^4)./4+(a5.*T.^5)./5+a6);
 else
     h = R*(b1.*T+(b2.*T.^2)./2+(b3.*T.^3)./3+(b4.*T.^4)./4+(b5.*T.^5)./5+b6);
 end

For specific heat calculation

% Calculation specific heat
function cp=specificheat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,global_med_temp)
%comparing with globaltemp
% If less than global it will take min coeffs
if T > global_med_temp
    cp = R.*(a1+a2.*T+a3.*T.^2+a4.*T.^3+a5.*T.^4);
 % If temp is more it will take max coeffs
else
    cp = R.*(b1+b2.*T+b3.*T.^2+b4.*T.^3+b5.*T.^4);
end

https://www.youtube.com/watch?v=xH2ArJK2tO8&t=11s