Project 1 - Parsing NASA thermodynamic data

AIM : To write a code in MATLAB to parse the NASA thermodynamic data file and then calculate thermodynamic properties of various gas species.

Objective:

1. To write a function that extracts the 14 coefficients and calculates the enthalpy, entropy and specific heats for all the species in the NASA thermodynamic data file "THERMO.dat"

The NASA polynomials have the form as shown below

Where R is the universal gas constant, T is the temperature.

2. To calculate the molecular weight of each species and display it in the command window.

3. To plot Cp, Enthalpy and Entropy for the local temperature range(low temperature:high temperature)specific for each species.

4. To save the plots as images with appropriate names and dump them in seperate folders for each species with suitable name automatically.

File Parsing:

File parsing refers to how information can be read from a particular text file or any type of file and how that information can be used in a computer program.

The term has slightly different meanings in different branches of linguistics and computer science. Traditional sentence parsing is often performed as a method of understanding the exact meaning of a sentence or word, sometimes with the aid of devices such as sentence diagrams. It usually emphasizes the importance of grammatical divisions such as subject and predicate.

Within computational linguistics the term is used to refer to the formal analysis by a computer of a sentence or other string of words into its constituents, resulting in a parse tree showing their syntactic relation to each other, which may also contain semantic and other information.

NASA Thermodynamic File:

THERMO.dat file is a data file containing 14 coefficients for 53 different elements given by NASA. These coefficients are utilized to calculate the specific heat,enthalpy and entropy for a particular element. 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 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 N20 NCO NH NH3 NNH NO NO2 O O2 OH

Enthalpy,Entropy and Specific Heat :

Enthalpy, entropy and specific heaat of various species values at different temperatures are required to analyze the thermodynamic process.

NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp,H,S. They have documented the coefficients that are required to evaluate these polynomials for 1000+ species.

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 coefficients 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 coefficients (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 coefficients (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).

Specific Heat:

The specific heat is the amount of heat per unit mass required to raise the temperature by one degree Celsius. The relationship between heat and temperature change is usually expressed in the form shown below

 Cp = R(a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4)

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.

 S  = R(a1 lnT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7)

Enthapy:

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 its value is determined entirely by the temperature, pressure, and composition of the system and not by its history.

 H = RT(a1 + a2 T /2 + a3 T^2 /3 + a4 T^3 /4 + a5 T^4 /5 + a6/T)

Main Code:

clear all
close all
clc

%inputs
R=8.314; %Gas Coefficient
%Reading data file
%open file and reading THERMO.dat file
f1=fopen('THERMO.DAT','r')
get=fgetl(f1);

%Reading Temperature from the file
temp= fgetl(f1);

%splitting
s=strsplit(temp);

%Converting the split string cells into numbers

global_low_temp = str2num(s{2});
global_medium_temp = str2num(s{3});
global_high_temp = str2num(s{4});

%Temperature Range
Temp=linspace(global_low_temp ,global_high_temp,1000);

%Skipping Comment lines
for i=1:3
    get=fgetl(f1);
end

%Reading the species
for i=1:53
    tline =fgetl(f1);
    
    A=strsplit(tline,' ');
    species=(A{1});

%Converting string to number
lowtemp =str2num(A{end-3});
midtemp =str2num(A{end-2});
hightemp =str2num(A{end-1});

tline1=fgetl(f1); %getting next line
a=strfind(tline1,'E'); %finding E

%Finding first 7 higher coefficients andsecond seven lower coefficients
%Finding Higher Coefficients
a1=tline1(1:a(1)+3);
a1=str2num(a1);

a2=tline1(a(1)+4:a(2)+3);
a2=str2num(a2);

a3=tline1(a(2)+4:a(3)+3);
a3=str2num(a3);

a4=tline1(a(3)+4:a(4)+3);
a4=str2num(a4);

a5=tline1(a(4)+4:a(5)+3);
a5=str2num(a5);

tline2 = fgetl(f1);
b =strfind(tline2,'E');
a6=tline2(1:b(1)+3);
a6=str2num(a6);

a7=tline2(b(1)+4:b(2)+3);
a7=str2num(a7);

%Finding lower coefficients

a8=tline2(b(2)+4:b(3)+3);
a8=str2num(a8);

a9=tline2(b(3)+4:b(4)+3);
a9=str2num(a9);


a10=tline2(b(4)+4:b(5)+3);
a10=str2num(a10);

tline3=fgetl(f1); %Getting to next line
c=strfind(tline3,'E');

a11=tline3(1:c(1)+3);
a11=str2num(a11);

a12=tline3(c(1)+4:c(2)+3);
a12=str2num(a12);

a13=tline3(c(2)+4:c(3)+3);
a13=str2num(a13);

a14=tline3(c(3)+4:c(4)+3);
a14=str2num(a14);

%Calculate Specific heat

Cp= Specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp);

%Calculate Entropy

S= Entropy_calc(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,Temp,R,global_medium_temp);

%Calculate Enthalpy

H= Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,Temp,R,global_medium_temp);

%Calculate Molecular Weight

Molecularweight = Molecular_weight(species);

%Writing molecular weight file

MW_file = fopen('Molecular Weight.txt','a+');
fprintf(MW_file,'Molecular weight of %s is %u \n',species,Molecularweight);
fclose(MW_file); %Closing the file

%Plotting Cp,H and S with Temperature
%Saving in a folder
cd('C:\Users\mohit\OneDrive\Desktop\project1\Species') 
mkdir(species);
cd(species);


%Plot1 Cp Vs T
figure(1)
plot(Temp,Cp,'linewidth',2,'color','b');
xlabel('Temperature(K)');
ylabel('Specific Heat(KJ)');
title(sprintf('Specific heat vs Temperature range of %s',species));
grid on
saveas(figure(1),'Specific Heat.jpg');


%Plot2 S Vs T
figure(2)
plot(Temp,S,'linewidth',2,'color','r');
xlabel('Temperature(K)');
ylabel('Entropy(KJ/mol-K');
title(sprintf('Entropy vs Temperature range of %s',species));
grid on
saveas(figure(2),'Entropy.jpg');

%Plot3 H Vs T
figure(3)
plot(Temp,H,'linewidth',2,'color','g');
xlabel('Temperature(K)');
ylabel('Enthalpy(KJ/mol)');
title(sprintf('Enthalpy vs Temperature range of %s',species));
grid on
saveas(figure(3),'Enthalpy.jpg');

cd('C:\Users\mohit\OneDrive\Desktop\project1') %Changing directory to original MATLAB path






end

Code Explanation:

• Clear all, Close all, clc commands are used to clear all the variables in workspace , close all the plots and figures and to clear the command window.
• Firstly, we will input the parameter R and we will open a file named 'THERMO.dat' and read this file using fopen command.
• Then, we read the lines in it using fgetl command and then we read the temperature file and we split the temperatures using strsplit command to convert it into cell array and then convert this string to number using str2num command to find low,high and mid temperatures.
• Using 'linspace' specified the temperature range from low to high.
• Using 'for and end' loops to skip the next three comment lines in the text file.
• Reading all the 53 species to plot specific heat(cp),enthalpy(H) and entropy(S) with temperature values and to fing molecular weights.
• Used 'for' loop to read 53 species.
• Used 'Strfind' to find the position of string.
• Extracted higher and lower temperature coefficients by reading line by line from the species.
• Out of which we will find first 7 coefficients which are for higher temperatures and the second seven coefficients which are for lower temperature.
• Calling specific heat, enthalpy, entropy function codes.
• calculating molecular weight from the molecular weight function.
• writing the molecular weight.txt file and fprintf all the molecular weights of the species and then close the txt file using fclose command.
• Change the directory using cd command and then creating new directory folder where we save the files using mkdir command and then change directory to this created folder.
• Plot Specific heat, entropy, enthalpy with the temperature for different species. Used 'Scanf' command inside the 'title'command to create different species title.
• Used 'Saveas' command to save all the plots of different species in the directory folder and change the directory to which our main code file runs in MATLAB.

Function to find Specific Heat :

% Specific heat calculation function

function Cp = Specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp)
%comparing with global temperature
%If less than global temperature minimum coefficients are taken
if Temp <= global_medium_temp
    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));
end

Explanation:

• We found Specific heat by creating specefic heat function and assigning its related formula coefficient from 14 coefficients, Gas constant, Temperature and mid temperature.
• Created a condtion using if loop and made the temperature less than or equal to mid temperature to specify lower temperature coefficient in the formula.
• Used else condition to specify coefficient of higher temperature in formula.
• End the loop to create a function to calculate Specific Heat.

Function to find 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
     S= R*(a1.*log(Temp)+ (a2.*Temp)+(a3*((Temp).^2)/2)+(a4*((Temp).^3)/3)+(a5*((Temp).^4)/4)+a7);
end

Explanation:

• We found Entropy by creating entropy function and assigning its related formula coefficient from 14 coefficients, Gas constant, Temperature and mid temperature.
• Created a condtion using if loop and made the temperature less than or equal to mid temperature to specify lower temperature coefficient in the formula.
• Used else condition to specify coefficient of higher temperature coefficient in formula.
• End the loop to create a function to calculate Entropy.

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

Explanation:

• We found Enthalpy by creating Enthalpy function and assigning its related formula coefficient from 14 coefficients, Gas constant, Temperature and mid temperature.
• Created a condtion using if loop and made the temperature less than or equal to mid temperature to specify lower temperature coefficient in the formula.
• Used else condition to specify coefficient of higher temperature coefficient in formula.
• End the loop to create a function to calculate Enthalpy.

Function to calculate Molecular Weight of different Species :

%Molecular Weight Calculation function
function Molecularweight =Molecular_weight(species)
% standard atomic name
atomic_name=['H','C','O','N','A'];
%atomic weight of Hydrogen,Carbon,Oxygen,Nitrogen,Argon

atomic_weight=[1.008 12.011 15.999 14.007 39.948];

Molecularweight=0 ; % Molecular Weight starting range

for i=1:length(species)  
    for j=1:length(atomic_name)
        if strcmp(species(i),atomic_name(j))
            Molecularweight = Molecularweight + atomic_weight(j);
            position = j;
        end
    end
    n= str2num(species(i));
    if n>1
        Molecularweight = Molecularweight + atomic_weight(position)*(n-1);
    end
end
fprintf('Molecularweight of %s :',species);
fprintf('%d',Molecularweight);
disp('  ');
end

Explanation:

• Firstly, created a Molecular_weight function and assigning species input to it.
• Defined the main elements of all 53 species in an array atomic_name and those are "H","C","O","N","Ar".
• Taking main elements array, created atomic weight array.
• Assigned molecular weight range starting from zero and then created a loop condition which runs through the length of species and another nested loop is created which runs which runs through the length of string atomic_name.
• Created a condition using if condition and used strcmp command to compare the species and atomic name of the strings to calculate molecular masses if by comaparision both strings are same, then assigned to its molecular weight.
• Used str2num of species and assigned its value to n. Ran a condition of if loop to find whether species contains more than one element.
• If more than one element i.e n>1, then incremented the atomic weight by multiplication of position and n-1 to atomic_weight array.
• Printed the Molecular Weight using fprintf command and displayed it using disp command.

Errors:

• Typo error in the above three, corrected them to fgetl,Mw_file,linewidth

• This error is faced due to invalid syntax, that is assigned the cell arrays by using (). Corrected this eror by assigning curly bracket in place of ().

Final Results:

Plots:

plots of O2 species :

Plots of N2 Species :

Plots of CO2 Species:

Screenshot of directory where different species are saved:

Outputs of Molecular weight of different species :

Conclusion:

Using MATLAB, the necessary coefficients required to evaluate various thermodynamic parameters for different species are extracted from NASA thermodynamic data file. The specific heat, enthalpy, entropy for all the species are calculated and plotted against the local temperature range. The molecular weight is calculated for all the species.

Drive Link for folders:

https://drive.google.com/drive/folders/1DMlOXS_roBycU1cPnBMbZvDbGO4rh7BU?usp=sharing

References:

https://in.mathworks.com/help/matlab/ref/fprintf.html

https://in.mathworks.com/help/matlab/ref/mkdir.html

https://en.wikipedia.org/wiki/Parsing

http://combustion.berkeley.edu/gri-mech/data/nasa_plnm.html