Project 1 - Parsing NASA thermodynamic data

Aim – Write a code in MATLAB for Find Cp, H and S from a thermo data file

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.
5. The main program that you will submit will just take the "THERMO.dat" file as an input and generate Cp, H and S plots.

NASA Thermodynamic File

File parsing is a string or symbol either in natural language computer language or data structure conferming to the rules of the former grammer.

THERMO.dat file is a data file that contain 14 coefficient for 53 different elements which file given by NASA. These coefficients are utilized to calculate the specific heat, enthalpy & entropy. Each element has its own range of local temperature.

List of 53 Elements contained in 'THERMO.dat' file:

AR C C2H C2H2 C2H3 C2H4

C2H5 C2H6 C3H7 C3H8 CH CH2

CH2(S) CH2CHO CH2CO CH2O CH2OH CH3

CH3CHO CH3O CH3OH CH4 CN CO

CO2 H H2 H2CN H2O H2O2

HCCO CCOH HCN HCNN HCNO HCO

HNCO HNO HO2 HOCN N N2

N2O NCO NH NH2 NH3 NNH

NO NO2 O O2 OH

MATHEMATICAL EQUATION:

For specific heat

`(C_p)/R = a_1 + a_2T + a_3T^2 + a_4 T^3 + a^5 T^4`

For enthalpy

`H/(RT) = a_1 + a_2T+ a_3T^2/3 +a_4T^3/4 +a_5T^4/5+a^6/T`

For entropy

`S/R = a_1 logT + a_2T + a_3T^2/2 + a^4T^3/3 +a^5T^4/4 +a_7`

where

T=local temprature

a1-a7=high temperature coefficient

a8-a14=low temprature coefficient

R=Gas constant=8.314

PROGRAM:

The program for file parsing

clear
close
clc

mkdir 'C:\Users\rajpu\Onedrive\Desktop\NASA'

%For Reading data file
    f1=fopen('thermo.dat','r');
    gt= fgetl(f1);

%For Finding global temprature
    global_temp=fgetl(f1);

%Splitting Coefficient
    g_temp=strsplit(global_temp);

%Converting string cells into number using str2num
    global_low_temp   = str2num(g_temp{2});
    global_medium_temp= str2num(g_temp{3});
    global_high_temp  = str2num(g_temp{4});

%Temprature Range
    Temp=linspace(global_low_temp,global_high_temp,1000);
% Gas Constant
    R=8.314;

% using loop for skipping 3 lines
    for i = 1:3
        gt=fgetl(f1);
    end

    for i = 1:53;
        gt=fgetl(f1); %Printing by skipping 1 line
        A=strsplit(gt); %Splitting
        Species_SP=(A{1}); %select one char


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

% %  Extracting higher and lower temperature coefficients
% 
% %  Extracting FIRST 7 coefficients i.e HIGH-temperature 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 SECOND 7 coefficients - LOW-temperature 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_calc(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,Temp,R,global_medium_temp);

    % Calculating molecular weight of species
    Molecularweight(i) = Molecular_weight(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)
        % Saving as a folder
 
           cur_dir=pwd;
           mkdir(Species_SP);
           cd(Species_SP);
 
        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 = Molecular_weight(Species_SP)
        f2 = fopen('molecular mass','w');
        fprintf('%s n %s n %f n','species','Molecular mass',M);
        fclose(f2)
end

Function Code For Specific Heat

% Specific heat calculation function
 
function Cp = Specificheat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,Medium_temp_range )          
% Comparing with global temperature
% If less than global temperature - minimum coefficients are taken

if  Temp <= Medium_temp_range 
    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));
    % Else maximum coefficients are taken
end

Function Code For Enthalpy

% Enthalpy Calculation function

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
 
       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
% Else maximum coefficients are taken 
           
        H = R*Temp.*(a1 + ((a2*Temp)/2) + ((a3*(Temp).^2)/3) + ((a4*(Temp).^3)/4) + ((a5*(Temp).^4)/5) +((a6./Temp)));
end

Function Code For Entropy

% Entropy Calculation function

function S = Entropy_calc(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
 
    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
 
% Else - maximum coefficients are taken
        S = R.*(a1.*log(Temp) + (a2.*Temp) + (a3*((Temp).^2)/2) + (a4*((Temp).^3)/3) + (a5*((Temp).^4)/4) + a7);
end

Function Code For Molecular Weight

% Molecular Weight Calculation function

function Molecularweight = Molecular_weight(Species_SP)
         % Standard - Atomic Name
         atomic_name = ['H','C','O','N','A'];             
         % Hydrogen, Carbon, Oxygen, Nitrogen, Argon
         
         fprintf('Molecular Weight of ');
         fprintf(Species_SP);
         fprintf(' : ');
         % Standard Atomic Number
         Atomic_Weight = [1.008 12.011 15.999 14.007 39.948];     
         % Starting Range of Molecule           
         Molecularweight = 0;                                   
        
         % Adding Atomic number to obtain 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 %s :',Species_SP);
            fprintf('%d',Molecularweight);
            disp(' ');      
end

Code Explanation:-

Clear close CLC command is used for clear screen clear work space and clear command window

mkdir command is used for making new folder

fopen command is used for read a file

fgetl command used for read a line from a data file

Strsplit used for split string or character vector

Str2num convert character array or string to numeric array

Space is used for dividing a line into a specific number of part

R is a gas constant

Now using for loop for skipping three line

Now using for loop for repeating a program 53 times for every coefficient of element

Strfind is used for find string within other string

Need to dividing the cooefficient and finding a number before E

Now writing a program to extract a coefficient from all lines coefficient or in the numbers of a1 a2 a3 .....

Creating a function for calculating for specific heat calculation, enthalpy calculation and Entropy calculation and Molecular weight calculation

Now creating a plot for a different -different condition for all the elements in of a data file and printing a molecular weight by using a command a fprintf command.

OUTPUT:

ScreenShot where all the folder are saved:

Screenshot of Molecular Weight

Attachment Of Plots for O2

Attachment Of Plots for N2

Attachment Of Plots for CO2