Project 1 - Parsing NASA thermodynamic data

Aim:

To write a MATLAB program to parse the NASA thermodynamic data file and calculate specific heat,enthalpy,entropy, molecular weight and plot the nature for various species.

Objectives:

• To write a function to extract 14 coeffs of various species from given thermodynamic data
• To write a function to extract enthalpy , entropy, specific heat using NASA polynomials of given species
• To write a function that plots specific heat, entropy, enthalpy of given species
• To write a function that calculates the molecular weight of given species
• To save the molecular weight data in seperate folders

Theory:

Parsing means breaking the datablock into smaller pieces so it can be easily interpreted.

• First give the input in the program for universal gas constant R = 8.314.
• The given data is formatted as the therm.dat file
• Use the command 'fopen' and the keyword 'r' to read the NASA Thermodynamic file
• Use keyword fgetl to read line by line
• Ignore the title line
• From next line extract the global max, min , mid temperatures
• Use for loop to ignore the next three comment lines
• Then use a for loop to extract the 53 species of data given in the file and calculate the required values
• Use fgetl to read the next line inside the loop and use strtok command to store the name of the element
• use str2num command to convert the string values that are read from next 3 lines to numbers that can be used for further calculations
• The low temperature values are stored in variables a1,a2,a3,a4,a5,a6,a7 and high temperature values are stored in variables a8,a9,a10,a11,a12,a13,a14
• Then use the references for the functions of specific heat, enthalpy, entropy, plots, molecular weight and write the molecular weights of the different elements to a file
• Then use fclose command to close the file

clear all
close all
clc

%Input
R = 8.314;

%Read the file
f1 = fopen('THERMO.dat','r');

%Read the first line of file
L = fgetl(f1);

%read and get global minimum,maximum and mid temperatures
L1 = fgetl(f1);
t = strsplit(L1);
global_min = str2num(t{2});
global_mid = str2num(t{3});
global_max = str2num(t{4});

%Define temperture range
temp = linspace(global_min,global_max,500);

%read and ignore the comment lines
for i = 1:3
    fgetl(f1);
end

for i = 1:53
    %read the first line to get species information
    l=fgetl(f1);
    % read the species name
    element_name = strtok(l);
    %read the lines for 7 high and 7 low temperature coeffs
    l1 = fgetl(f1);
    l2 = fgetl(f1);
    l3 = fgetl(f1);
    %Find E to extract the coefficients
    a= strfind(l1,'E');
    b= strfind(l2,'E');
    c= strfind(l3,'E');

    %Extract high temp coeffs
    a1 = str2num(l1(1:a(1)+3));
    a2 = str2num(l1(a(1)+3:a(2)+3));
    a3 = str2num(l1(a(2)+4:a(3)+3));
    a4 = str2num(l1(a(3)+4:a(4)+3));
    a5 = str2num(l1(a(4)+4:a(5)+3));
    a6 = str2num(l2(1:b(1)+3));
    a7 = str2num(l2(b(1)+4:b(2)+3));

    %Extract low temp coeffs
    a8 = str2num(l2(b(2)+4:b(3)+3));
    a9 = str2num(l2(b(3)+4:b(4)+3));
    a10 = str2num(l2(b(4)+4:b(5)+3));
    a11 = str2num(l3(1:c(1)+3));
    a12 = str2num(l3(c(1)+4:c(2)+3));
    a13 = str2num(l3(c(2)+4:c(3)+3));
    a14 = str2num(l3(c(3)+4:c(4)+3));

    %Calculate specific heat
    c_p = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,temp,global_min,global_max,global_mid);

    %Calculate enthalpy
    H = enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,temp,global_min,global_max,global_mid);

    %Calculate entropy
    S = entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,temp,global_min,global_max,global_mid);

    %Plotting graphs
    plots(c_p,H,S,temp,element_name)

    %Calculate Molecular weight
    mol_wt(i) = molecular_weight(element_name);
    disp(mol_wt(i));
    Molecular_Weight = strcat(mol_wt(i),element_name);

    %Function to write molecular weight to a file
    weight = mol_wt(i);
    writing_molwt(element_name,weight);
end
fclose(f1);

Function to calculate the specific heat:

• This functions use the temperature coefficients, gas constant and global temperature values
• Then use if else statement to check if the temperature value is less global_mid and the appropriate temperature coefficients to calculate the c_p value

%Function to calculate specific heat Cp

function c_p = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,temp,global_min,global_max,global_mid)
if temp<=global_mid
    c_p = R.*(a1 + a2.*temp + a3.*(temp.^2) + a4.*(temp.^3) + a5.*(temp.^4));
else
    c_p = R*(a8 + a9.*temp + a10.*(temp.^2) + a11.*(temp.^3) + a12.*(temp.^4));
end

Function to calculate the enthalpy:

• This function use the temperature coefficients, gas constant and global temperature values
• Then use if else statement to check if the temperature value is less global_mid and the appropriate temperature coefficients to calculate the H value

% Function to calculate the enthalpy 
function H = enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,temp,global_min,global_max,global_mid)
if temp <= global_mid
    H = R.*(a1 + a2.*(temp./2) + a3.*((temp.^2)./3) + a4.*((temp.^3)./4) + a5.*((temp.^4)./5) + a6./temp);
else
    H = R.*(a8 + a9.*(temp./2) + a10.*((temp.^2)./3) + a11.*((temp.^3)./4) + a12.*((temp.^4)./5) + a13./temp);
end 

Function to calculate the entropy:

• This function use the temperature coefficients, gas constant and global temperature values
• Then use if else statement to check if the temperature value is less global_mid and the appropriate temperature coefficients to calculate the S value

%Function to calculate entropy
function S = entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,temp,global_min,global_max,global_mid)
if temp <= global_mid
    S = R * (a1.*log(temp) + a2.*temp + a3.*((temp.^2)./2) + a4.*((temp.^3)./3) + a5.*((temp.^4)./4) + a7);
else
    S = R * (a8.*log(temp) + a9.*temp + a10.*((temp.^2)./2) + a11.*((temp.^3)./3) + a12.*((temp.^4)./4) + a14);
end

Function to draw the plots

• This function is used to draw plots temp vs specific heat,temp vs enthalpy, temp vs entropy
• Using mkdir command new folder is created in current folder with the same name as the element
• Then it is saved respectively in a folder
• Then cd command is used to go back to original folder

% Function to plot graph for elements
function plots(c_p,H,S,temp,element_name)
cd(pwd);
dir_name = element_name
mkdir(dir_name)

%plot temp vs specific heat
figure(1)
plot(temp,c_p,'LineWidth',3,'Color','r')
title(sprintf('Specific Heat vs Temperature %s:',element_name));
xlabel('Temp');
ylabel('Specific Heat');

%plot temp vs enthalpy
figure(2)
plot(temp,H,'LineWidth',3,'Color','g');
title(sprintf('Enthalpy vs Temperature %s:',element_name));
xlabel('Temp');
ylabel('Enthalpy');

%plot temp vs entropy
figure(3)
plot(temp,S,'LineWidth',3,'Color','b');
title(sprintf('Entropy vs Temperature %s:',element_name));
xlabel('Temp');
ylabel('Entropy');
cd(dir_name);

saveas(1,'cpvstemp.png');
saveas(2,'Hvstemp.png');
saveas(3,'Svstemp.png');

cd ..;

end

Function to molecular_weight :

• This function is used to calculate the molecular weight from the atomic weight
• The atomic weight is given in as array and strcmp command between element name and atomic weight of each element

%Function to calculate molecular weight
function [w] = molecular_weight(element_name)
s = {'H' 'C' 'N' 'O' 'AR'};

%Atomic weight of the elements
atomic_weight = [1 12 14 16 40];
w = 0;
for i =1:length(element_name)
    for j = 1:length(s)
        if strcmp(element_name(i),s(j))
            w = w + atomic_weight(j);
            %storing the position of the compared and matched element
            position = j;
        end
    end
    z = str2num(element_name(i));
    if z > 1
        w = w + atomic_weight(position)*(z-1);
    end
end
end

Function to write the molecular weight to file:

• This function is used to write the molecular weight of the element to the seperate file using append command
• The a+ command is used to append the molecular weight values in same file one after other

%Function to write molecular weights to a file
function writing_molwt(element_name,weight)

dir_name = element_name;
mol_file = fopen('mol_weight.dat','a+');
fprintf(mol_file,'Molecular weight of %s is %f n',dir_name,weight);
fclose(mol_file);
end

Conclusion:

Thus the NASA data is parsed successfully and the thermodynamic properties, molecular weight are calculated and stored in the respective folders successfully.Please find the link fro graphs and molecular weight below

https://drive.google.com/drive/folders/1ADMB-c_bqyUddI_oH0ruWZ8kd0hr2tNL?usp=sharing