Project 1 - Parsing NASA thermodynamic data

AIM:

Write a MATLAB code to calculate specific heat at constant pressure (Cp), Enthalpy (H) and Entropy (S) of 53 species of gases given in the data along with thier molecular weight.

GOVERNING EQUATIONS:

1. Cp = (a1 + a2*T + a3*T^2 + a4*T^3 + a5*T^4)*R
2. H = (a1 + a2*T/2 + a3*T^2/3 + a4*T^3/4 + a5*T^4/5 + a6/T)*R*T
3. S = (a1*ln(T) + a2*T + a3*T^2/2 + a4*T^3/3 + a5*T^4/4 + a7)*R

Cp = Specific heat at constant pressure (JK-1kg-1)

H = Enthalpy (J)

S = Entropy (JK-1)

R = Universal gas constant (Jmol-1K-1)

T = Local temperature (K)

OBJECTIVE OF PROJECT:

1. Extract coefficients and the local temperature specific to different species from NASA thermodynamic datasheet on 53 species of gases
2. Calculate specific heat (Cp), enthalpy (H) and entropy (S) for local temperature specific species
3. Plot graphs of specific heat (Cp), enthalpy (H) and entropy (S) against local temperature
4. Calculate molecular weight and display it in the command output
5. Save the plot as images with appropriate names and dump them in separate folders for each species with suitable name generated automatically

BODY OF THE PROJECT:

A MATLAB CODE

A function used to calculate Specific heat, Enthalpy and Entropy

%Function to calculate the values of Cp, H and S
function [Cp, H, S] = func_cal(s11,s12,s13,s14,s15,s21,s22,s23,s24,s25,s31,s32,s33,s34,TR,Local_mid_temp)
        
        %Cp is specific heat at constant pressure (JK-1kg-1)
        %H is the enthalpy (J)
        %S is the entropy (Jk-1)
        R = 8.314; %SI unit is Jmol^-1K^-1
        
        for k = 1:length(TR)
            
                if(TR(k)<=Local_mid_temp)
                   
                    Cp = (s23+(TR).*s24+(TR).^2.*s25+(TR).^3.*s31+(TR).^4.*s32).*R;
                    H = (s23+(s24.*TR)/2+(s25.*(TR).^2)/3+(s31.*(TR).^3)/4+(s32.*(TR).^4)/5+s33./TR).*R.*TR;
                    S = (s23.*log(TR)+s24.*TR+s25.*(TR).^2/2+s31.*(TR).^3/3+s32.*(TR).^4/4+s34).*R;
                
                else
                    
                    Cp = (s11+(TR).*s12+(TR).^2.*s13+(TR).^3.*s14+(TR).^4.*s15).*R;
                    H = (s11+(s12.*TR)/2+(s13.*(TR).^2)/3+(s14.*(TR).^3)/4+(s15.*(TR).^4)/5+s21./TR).*R.*TR;
                    S = (s11.*log(TR)+s12.*TR+s13.*(TR).^2/2+s14.*(TR).^3/3+s15.*(TR).^4/4+s22).*R;
                    
                end
                 
            
        end

end

A function to calculate molecular weight of the different species

%function to calculate the molecular weight of the species
function mol_wt = molecular_weight(species)
   S_atomic = {'H' 'C' 'N' 'O' 'AR'};
   at_wt = [1 12 14 16 40];
   weight = 0;
   
   for q = 1:length(species)
      for p = 1:length(S_atomic)
        if strcmp(species(q), S_atomic(p))==1
            weight = weight+at_wt(p);
            position = p;
        end
      end
      n = str2double(species(q));
     
      if n>1
         weight = weight + at_wt(position)*(n-1);
      end
   end
   

end

MAIN PROGRAM

%Parsing NASA thermodynamic data
clear all
close all
clc

%To use the read mode and read the title
thermal_data1 = fopen('thermal_data.dat', 'r');
file_name = fgetl(thermal_data1);
first_line = fgetl(thermal_data1);

%To read the global temperatures available
A = strsplit(first_line, ' ');
Global_low_temp = str2double(A{2});
Global_mid_temp = str2double(A{3});
Global_high_temp = str2double(A{4});

%To read the next three lines which are under comented sections
for i = 1:3
    
    C = fgetl(thermal_data1);
    
end

%To extract the coefficients of 53 individual species
for j = 1:53
    
    local_temp = fgetl(thermal_data1);
    T1 = strsplit(local_temp, ' ');
    species = (T1{1});
    temp = strfind(local_temp, 'G   ');
    
    Local_low_temp = str2double(local_temp(temp+1:temp+7));
    Local_mid_temp = str2double(local_temp(temp+21:temp+27));
    Local_high_temp = str2double(local_temp(temp+11:temp+16));
    TR = linspace(Local_low_temp, Local_high_temp, 500);
    
    s1line = fgetl(thermal_data1);
    b = strfind(s1line, 'E');
    s11 = str2double(s1line(1:b(1)+3));
    s12 = str2double(s1line(b(1)+4:b(2)+3));
    s13 = str2double(s1line(b(2)+5:b(3)+3));
    s14 = str2double(s1line(b(3)+4:b(4)+3));
    s15 = str2double(s1line(b(4)+5:b(5)+3));
    
    s2line = fgetl(thermal_data1);
    b2 = strfind(s2line, 'E');
    s21 = str2double(s2line(1:b2(1)+3));
    s22 = str2double(s2line(b2(1)+4:b(2)+3));
    s23 = str2double(s2line(b2(2)+5:b(3)+3));
    s24 = str2double(s2line(b2(3)+4:b(4)+3));
    s25 = str2double(s2line(b2(4)+5:b(5)+3));
    
    s3line = fgetl(thermal_data1);
    b3 = strfind(s3line, 'E');
    s31 = str2double(s3line(1:b3(1)+3));
    s32 = str2double(s3line(b3(1)+4:b(2)+3));
    s33 = str2double(s3line(b3(2)+5:b(3)+3));
    s34 = str2double(s3line(b3(3)+4:b(4)+3));
    
    
   [Cp, H, S]  = func_cal(s11,s12,s13,s14,s15,s21,s22,s23,s24,s25,s31,s32,s33,s34,TR,Local_mid_temp);
         
   %To create a folder with species name
   folder_name = sprintf('[%d]    [%s]', j, species);
   current_dir = pwd;
   
   %All the function program and thermodynamic data file are stored in a folder, that specify current directory
   mkdir(folder_name);
   cd(folder_name);
   
   %To plot the graphs
   figure(1)
   plot(TR, Cp, 'linewidth', 1)
   grid on
   xlabel('Local Temperature [K]')
   ylabel('Specific Heat with constant Pressure [JK-1mol-1]')
   title('Specific Heat (Cp) vs Temperature')
   IMAGE1 = sprintf("Cp vs TR of [%s].png", species);
   saveas(gca, IMAGE1)
   
   figure(2)
   plot(TR, H, 'linewidth', 1)
   grid on
   xlabel('Local Temperature [K]')
   ylabel('Enthalpy [J]')
   title('Enthalpy [H] vs Local Temperature')
   IMAGE2 = sprintf("H vs TR of [%s].png", species);
   saveas(gca, IMAGE2)
   
   figure(3)
   plot(TR, S, 'linewidth', 1)
   grid on 
   xlabel('Local Temperature[K]')
   ylabel('Entropy')
   title('Entropy [S] vs Local Temperature')
   IMAGE3 = sprintf("S vs TR of [%s].png", species);
   saveas(gca, IMAGE3)
   
   cd(current_dir);
   
   mol_wt = molecular_weight(species); %function call of molecular weight calculation
    
   fprintf('Molecular weight of %s is %d.', species, mol_wt)
   fprintf('\n \n')
   
end

PROCEDURE:

1) First the new editor window is opened, the commands "clear all", "close all" and "clc" are used to clear the command window and the previously done programs

2) The aim of the program is commented using the "%" symbol on the top of the editor window

3) The program written is saved usin the extension '.m'

4) Using command "fopen" and in read mode the data to be used is read.

The syntax of fopen command is

fileID = fopen(filename, permission);

This command opens file with the type of access specified by permission

5) Using command "fgetl" the next line is read.

The syntax is

tline = fgetl(fileID);

6) "str2double" or "str2num" is used to convert the text in string into numbers or numerical values with precision.

The syntax is

X = str2double(str) or X = str2num(str)

7) "strsplit" is used to split the string at delimiters specified

The syntax is

c = strsplit(str, delimiter)

8) Using a "for-loops" 53 species data is read using the above commands within the for loop.

9) The calculation of specific heat, enthalpy and entropy is done in a different editor by creating a function which is called within the loop. The formulae given above are used to calculate for each of the species using their specific local temperatures. For temperatures above middle temperature the first 7 coefficients read from the datasheet are used, while for the remaing set of temperatures the next 7 coefficients read from the datasheet are used.

10) Now the values are plotted against the local temperatures and stored in the same directory in different folders using the command "cd()"

11) In a new editor another function is created to calculate the molecular weight of the individual species. Here initially each atom's elemental formulae is stored under the name "S_atomic" and their corresponding atomic weight is stored in "at_wt" variable. Using "if strcmp( )" loop is used to compare the "species formulae" and the "S_atomic" variable in order to calculate the value of molecular weight.

12) Using command "fprintf" the output value of the molecular weight is seen in the command window.

ERRORS and CORRECTION:

1) The molecular weight of elements like O2 is incorrect.

Corrected it by putting

n = str2double(species(q));

if n>1
mol_wt = mol_wt + at_wt(position)*(n-1);
end

out of the for loop having range "1:length(S_atomic)" as per my program code but within the for loop containing the range "1:length(species)".

2) I got an error of giving the wrong indexing value which comes like this

correction: This error appeared because "n = str2double(species(p));" contained p as it's indexing variable. p ranges from 1:length(S_atomic) which contains less value than species variable. "q" is the indexing variable which contains 53 which is same as the species. Hence the error gets corrected.

RESULT:

PLOTS FOR O2:

PLOTS FOR N2:

PLOTS FOR CO2:

MOLECULAR WEIGHT: