Parsing NASA Thermodynamic Data File and Obtaining Different Plots From it.

Parsing NASA Thermodynamic Data File

Aim:- To open 'Thermo.dat' file and use this data for plotting different graphs.

Theory:-

The NASA thermodynamic data file is parsed to extract useful information from the file, like species coefficients, temperature ranges. The data extracted is used to calculate the Specific Heat, Enthalpy and Entropy of a species in its temperature range and respective graphs are plotted.

The NASA polynomials have the form:
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

where a1, a2, a3, a4, a5, a6, and a7 are the numerical coeffcients supplied in NASA thermodynamic files. The first 7 numbers starting on the second line of each
species entry (five of the second line and the first two of the third line) are the seven coeffcients (a1 through a7, respectively) for the high-temperature range
(above 1000 K, the upper boundary is specified on the first line of the species entry). The following seven numbers are the coeffcients (a1 through a7,
respectively) for the low-temperature range (below 1000 K, the lower boundary is specified on the first line of the species entry).

Objective:-

To write a program which extract all the 14 co-efficients and calculates the enthalpy,entropy,and specific heat for all the elements in 'Thermo.dat' file.

Link for 'Thermo.dat' File:- https://drive.google.com/file/d/1AYMhoznYslGeY_yt29P7e-IMUjHbbtaO/view?usp=sharing

Code:-

factor function code

function [cp,h,s]=factor(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,T)

R=8.314;
for j=1:length(T)
    if T(j)>300
        
        cp=(a1+a2.*T+a3.*T.^2+a4.*T.^3+a5.*T.^4).*R;
        h=(a1+(a2.*T)./2+(a3.*T.^2)./3+(a4.*T.^3)./4+(a5.*T.^4)./5+(a6./T)).*R.*T;
        s=(a1.*log(T)+a2.*T+(a3.*T.^2)./2+(a4.*T.^3)./3+(a5.*T.^4)./4+a7).*R;
        
    else
        cp=(a8+a9.*T+a10.*T.^2+a11.*T.^3+a12.*T.^4).*R;
        h=(a8+(a9.*T)./2+(a10.*T.^2)./3+(a11.*T.^3)./4+(a12.*T.^4)./5+(a13./T)).*R.*T;
        s=(a8.*log(T)+a9.*T+(a10.*T.^2)./2+(a11.*T.^3)./3+(a12.*T.^4)./4+a14).*R;
        
    end
end
end


         

mole Function code

function [molecular_weight]=mole(element_name)

ELEMENT=upper(element_name);
constituents={'H','C','N','O','A'};
atomic_weight=[1 12 14 16 40];

molecular_weight=0;
for i=1:length(ELEMENT)
    for j=1:length(constituents)
        if strcmp(ELEMENT(i),constituents(j))
            molecular_weight=molecular_weight+atomic_weight(j);
        end
    end
    
    n=str2double(ELEMENT(i));
    
    if n>1
        molecular_weight=molecular_weight+atomic_weight(j)*(n-1);
    end
end

fprintf('molecular weight is:-%d',element_name)
fprintf('%d',molecular_weight)
fprintf('')
fprintf('amu')
fprintf('n')
    end

Main Program

clear all 
close all
clc

%Opening thermodynamic data file
F=fopen('THERMO.dat','r');
file_name=fgetl(F);
% Temperatures from data file
low_temp=300;
mid_temp=1000;
high_temp=5000;

for i=(1:3)
    comments=fgetl(F);
end
for i=1:53
    
    % Element Name
    element=fgetl(F);
    element_name=element(1:6);
    
    %Coefficients of first part
    line1=fgetl(F);
    l1=strfind(line1,'E');
    a1=str2double(line1(1:(l1(1)+3)));
    a2=str2double(line1((l1(1)+4):(l1(2)+3)));
    a3=str2double(line1((l1(2)+4):(l1(3)+3)));
    a4=str2double(line1((l1(3)+4):(l1(4)+3)));
    a5=str2double(line1((l1(4)+4):(l1(5)+3)));
    
    %Coefficients of second part
    line2=fgetl(F);
    l2=strfind(line2,'E');
    a6=str2double(line2(1:(l2(1)+3)));
    a7=str2double(line2((l2(1)+4):(l2(2)+3)));
    a8=str2double(line1((l2(2)+4):(l2(3)+3)));
    a9=str2double(line1((l2(3)+4):(l2(4)+3)));
    a10=str2double(line1((l2(4)+4):(l2(5)+3)));
    
    %Coefficients of third part
    line3=fgetl(F);
    l3=strfind(line3,'E');
    a11=str2double(line3(1:(l3(1)+3)));
    a12=str2double(line3((l3(1)+4):(l3(2)+3)));
    a13=str2double(line3((l3(2)+4):(l3(3)+3)));
    a14=str2double(line3((l3(3)+4):(l3(4)+3)));
   
    % Creating a Folder for plots
    mkdir(['E:Skill_lync',element_name]);
    
    T=linspace(low_temp,high_temp,500);
    
    % Thermodynamic Equation Function
    [cp,h,s]=factor(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,T);
    
    % Plot of Temperature Vs Entropy
    figure(1)
    plot(T,s,'linewidth',3,'color','c')
    title(element)
    xlabel('Temperature')
    ylabel('Entropy')
    
    % Plot of Temperature Vs Specific Heat
    figure(2)
    plot(T,cp,'linewidth',3,'color','b')
    title(element)
    xlabel('Temperature')
    ylabel('Specific Heat')
    
    % Plot of Temperature Vs Enthalpy
    figure(3)
    plot(T,h,'linewidth',3,'color','b')
    title(element)
    xlabel('Temperature')
    ylabel('Enthalpy')
    
   
    % Saving Plots
    saveas(1,'Entropy.png')
    saveas(2,'SpecificHeat.png')
    saveas(3,'Enthalpy.png')
    molecular_weight=mole(element_name);
end
fclose(F)

   
    

Explaination of Code:-

First the Thermo.dat file is opened

All the value are taken from the file and assigned to different variables

low temperature,mid temperature, high temperature are mentioned

specific heat,enthalpy,entropy are calculated using formula

different directory are made for each element

molecular weight, mole function is assigned for that

14 different coefficient are assigned with values

graphs are plotted of each element and stored in their respective folders.

temperature Vs entropy,temperature Vs specific heat,temperature Vs enthalpy are plotted and saved in png format.

Output:-

O2:-

N2:-

Folders formed:-

Error Occured:-

Due to this error it was unable to get all the graphs

Error can be rectified by using higher graphic card in you PC or Laptop

Conclusion:-

It can be concluded that enthalpy,entropy,specific heat and temperature are proportional to each other. Molecular weight changes with atomic weight.