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 coefficient that is required to evaluate these polynomials for 1000+ species.

Link to the File is given below -->

https://drive.google.com/file/d/1AikCHDlz9s_fu8whE81qee_GD4AyYpCP/view

The Function codes for -->

1. Entropy -->

Entropy, the measure of a system’s thermal energy per unit temperature that is unavailable for doing useful work. Because work is obtained from ordered molecular motion, the amount of entropy is also a measure of the molecular disorder, or randomness, of a system. The concept of entropy provides deep insight into the direction of spontaneous change for many everyday phenomena.

%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

Reference -->

https://in.mathworks.com/help/images/ref/entropy.html

2. Enthalpy -->

Enthalpy, the sum of the internal energy and the product of the pressure and volume of a thermodynamic system. Enthalpy is an energy-like property or state function—it has the dimensions of energy (and is thus measured in units of joules), and its value is determined entirely by the temperature, pressure, and composition of the system and not by its history. In symbols, the enthalpy, H, equals the sum of the internal energy, E, and the product of the pressure, P, and volume, V, of the system: H = E + PV.

%Calculation of enthalpy
function h = enthalpy(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

Reference -->

https://in.mathworks.com/matlabcentral/fileexchange/30026-enthalpy-calculator-gui

3. Specific Heat -->

%Calculation of specific heat
function cp = specificheat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,global_med_temp)
% Comparing with global temp
% 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

Reference -->

1. https://in.mathworks.com/help/physmod/simscape/ref/thermodynamicpropertiessensorg.html

2.

https://in.mathworks.com/help/pde/ug/pde.thermalmaterialassignment-properties.html

4. Molecular Weight -->

%Function to calculate molecular weight
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(i));
        if n>1
            MW = MW + atomic_weight(position)*(n-1);
        end
    end
    fprintf(\'%d\',MW);
    disp(\' \');
    end

Reference -->

1.

https://in.mathworks.com/matlabcentral/answers/106320-make-a-loop-to-calculate-molecular-weight

2.

https://in.mathworks.com/matlabcentral/fileexchange/26139-molecular-weight-calculator-function

The Main code for the main program that you will submit will just take the \"THERMO.dat\" file as an input and generate Cp, H and S plots for all species that are included in the file is,

clear all
close all
clc

%NASA\'s thermodynamic data

f1 = fopen(\'THERMO.dat\',\'r\');
first_line = fgetl(f1);

% Obtaining global max, min, and medium temperatures

global_temp = fgetl(f1);
temp_range = strsplit(global_temp,\' \');
global_min_temp = str2num(temp_range{2});
global_med_temp = str2num(temp_range{3});
global_max_temp = str2num(temp_range{4});

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

%skipping the next three comments in the file
for i=3:5
  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 lower 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);

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

Explanation to the code:

• The \'thermo.dat\' file is saved in the same directory were the MATLAB files are saved.
• The functions such as \'enthalpy’, \'entropy’, \'specific heat\' and ’molecular weight\' are defined using the equations stated.
• In the main code, necessary inputs such as universal gas constant value are assigned to variable R and the \'fopen\' command is used to open the \'thermo.dat\' file. Header information is extracted using \'fgetl\' command.
• Global temperature range is stored in the form of numbers using\'str2num\' command.
• For loop is used to read the data for all the 53 species one at a time. The required data from line0 to line3. Line0 describes the name of species and local temperature range for the species. Line1 to line3 contains 14 coefficients for the species out of which first 7 are high-temperature coefficients and the next 7 are low-temperature coefficients.
• These coefficients are extracted and converted from string data to numbers using \'str2double\' command.
• Then the local temperature range is assigned to variable t for the plots.
• The defined functions are called to calculate the corresponding properties of the species.

Function code for Plotting the data is,

%Create plotting 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(\'Enthalpy Vs Temperature for %s\',name))
xlabel(\'TEMPERATURE\')
ylabel(\'ENTHALPY\')

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\')
    
end

1. All the species are made of elements which are Hydrogen, Carbon, nitrogen, oxygen, Ar. So, for these their molecular weights had been taken as input.

2. Starting with for loop one for species name and other for elements, strcmp command is used to compare the elements with species if both matches then their molecular weight will be calculated alternately with every matching element with to species.

3. Finding a string in species name and converting into a number then taking their position for extracting their weight, if \'n\' it will be greater than one then it is used in a further step to calculate the molecular weight. this \'n\' here represents the number of molecules of an element in species.

The following results were obtained -->

Molecular Weight of O : 16
Molecular Weight of O2 : 32
Molecular Weight of H : 1
Molecular Weight of H2 : 2
Molecular Weight of OH : 17
Molecular Weight of H2O : 18
Molecular Weight of HO2 : 33
Molecular Weight of H2O2 : 34
Molecular Weight of C : 12
Molecular Weight of CH : 13
Molecular Weight of CH2 : 14
Molecular Weight of CH2(S) : 14
Molecular Weight of CH3 : 15
Molecular Weight of CH4 : 16
Molecular Weight of CO : 28
Molecular Weight of CO2 : 44
Molecular Weight of HCO : 29
Molecular Weight of CH2O : 30
Molecular Weight of CH2OH : 31
Molecular Weight of CH3O : 31
Molecular Weight of CH3OH : 32
Molecular Weight of C2H : 25
Molecular Weight of C2H2 : 26
Molecular Weight of C2H3 : 27
Molecular Weight of C2H4 : 28
Molecular Weight of C2H5 : 29
Molecular Weight of C2H6 : 30
Molecular Weight of CH2CO : 42
Molecular Weight of HCCO : 41
Molecular Weight of HCCOH : 42
Molecular Weight of H2CN : 28
Molecular Weight of HCN : 27
Molecular Weight of HNO : 31
Molecular Weight of N : 14
Molecular Weight of NNH : 29
Molecular Weight of N2O : 44
Molecular Weight of NH : 15
Molecular Weight of NH2 : 16
Molecular Weight of NH3 : 17
Molecular Weight of NO : 30
Molecular Weight of NO2 : 46
Molecular Weight of HCNO : 43
Molecular Weight of HOCN : 43
Molecular Weight of HNCO : 43
Molecular Weight of NCO : 42
Molecular Weight of CN : 26
Molecular Weight of HCNN : 41
Molecular Weight of N2 : 28
Molecular Weight of AR : 40
Molecular Weight of C3H8 : 44
Molecular Weight of C3H7 : 43
Molecular Weight of CH3CHO : 44
Molecular Weight of CH2CHO : 43
>>

Plots for O2 -->

Plots for N2 -->

Plots for CO2 -->

The error that was occurred to me during the programming were after the programming was complete and ran.

I noticed that the \'data \' file was not mentioned in the Directory of the program.

Which gave me this error