Challenge

Project 1 - Parsing NASA thermodynamic data

KKarun Das19 Aug 2022

Objective-

1) 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.

2) To calculate the molecular weight of each species

3) To plot the Cp, Enthalpy and Entropy for the local temperature range

Report

Equations used

Cp/R = a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4

H/RT = a1 + a2 T /2 + a3 T^2 /3 + a4 T^3 /4 + a5 T^4 /5 + a6/T

S/R = a1 lnT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7

Main program

% Main program
% Clearing workspace, plots, figures and command window
clear all
close all
clc

%  Reading Data file - NASA Thermodynamic
% Opening the data file using fopen
f1=fopen('THERMO.dat','r')

% Read line from file
gt=fgetl(f1);

% Assigning temperature values from the data file and assigning it to
% global_temp by reading the header tags

global_temp=fgetl(f1);

% Spliting
t=strsplit(global_temp);

% Converting string cells to numbers using str2num,
global_low_temp=str2num(t{2});
global_medium_temp=str2num(t{3});
global_high_temp=str2num(t{4});

% Temperature range & Gas Constant
% Temperature Range assumed using linspace command
Temp=linspace(global_low_temp,global_high_temp,1000);

% Universal Gas Constant J/(mol-K)
R = 8.314;   
    
% Skipping three comment lines in the file

for i = 1:3
    gt=fgetl(f1);
end

 for i = 1:53
gt=fgetl(f1);
A=strsplit(gt,' ');
Species_SP=(A{1})


line_1=fgetl(f1)  ;        
a=strfind(line_1,'E');

% Extracting higher and lower temperature coefficients

% Extracting FIRST 7 coefficients 

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 last 7 coefficients

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 calculation
    Cp = specificheat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp);
 
    % Enthalpy calculation
    H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,Temp,R,global_medium_temp);
 
    % Entropy calculation
    S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,Temp,R,global_medium_temp);

    % Calculating molecular weight of species
    Molecularweight(i) = molecularweight(Species_SP);
    
    % Read & Write Molecular list
     Molecular_list = fopen('Molecular weight.txt','a+');
         fprintf( Molecular_list,'Molecular weight of %s is %d \n',Species_SP,Molecularweight);
         fclose(Molecular_list);

     % plotting(CP,H,S,Temp,Species_SP
 
           cur_dir=pwd;
           mkdir(Species_SP);
           cd(Species_SP);
 
        % Plot 1, Specific heat (CP) vs Temperature Range (Temp)
 
        figure(1)
        plot(Temp,Cp,'linewidth',2,'color','r');
        xlabel('Temperature  [K]');
        ylabel('Specific Heat [kJ]');
        title(sprintf('Varying Specific Heat with Temperature Range of %s',Species_SP));
        grid on
        saveas(figure(1),'Specific Heat.jpg');
 
        % Plot 2, Enthalpy (H) vs Temperature Range (Temp)
 
        figure(2)
        plot(Temp,H,'linewidth',2,'color','g');
        xlabel('Temperature [K]');
        ylabel('Enthalpy in [kJ/(mol)]');
        title(sprintf('Varying Enthalpy with Temperature Range of %s',Species_SP));
        grid on
        saveas(figure(2),'Ethalpy.jpg');
 
        % Plot 3, Entropy (S) vs Temperature Range (Temp)
 
        figure(3)
        plot(Temp,S,'linewidth',2,'color','b');
        xlabel('Temperature [K]');
        ylabel('Entropy in [kJ/(mol-K)]');
        title(sprintf('Varying Entropy with Temperature Range of %s',Species_SP));
        grid on
        saveas(figure(3),'Entropy.jpg');
 
        cd(cur_dir);            
 
        % Molecular weight
        
        M = Molecularweight(Species_SP)
        f2 = fopen('molecular mass','w');
        fprintf('%s n %s n %f n','species','Molecular mass',M)
        fclose(f2)
end

Function definitions used

Specific heat

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

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 temperature
     % Minimum coefficients are taken
     % Else maximum 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

Entropy

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

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 temperature
       % Minimum coefficients are taken
       % Else maximum 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

Enthalpy

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

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 temperature
     % Minimum coefficients are taken
     % Else maximum 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

Molecular weight

function Molecularweight = molecularweight (Species_SP)
         
 % Atomic Name
   
 atomic_name = ['H','C','O','N','A'];             
     
 fprintf('Molecular Weight of ');
 fprintf(Species_SP);
 fprintf(' : ');
 
 % Standard Atomic Number
  Atomic_Weight = [1.008 12.011 15.999 14.007 39.948];                
  Molecularweight = 0;                                   
        
 
        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 %s :',Species_SP);
            fprintf('%d',Molecularweight);
            disp(' ');      
end

Result