Project 1 - Parsing NASA thermodynamic data

AIM : To Write a function that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the data file

Objective :

1. Calculate the molecular weight of each species and display it in the command window.
2. Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species.
3. Save the plots as images with appropriate names and dump them in separate folders for each species with suitable name. All these should be generated automatically.

File parsing : File parsing means the analysis of the required data file in a computer [for Ex THERMO.dat file] to get the Specific values which is required. In matlab we have a tool known as interactive data importing tool But in this project we will be coding a function to implement the file parsing as per the requirement .

NASA THERMO.dat file : The NASA THERMO.dat file will be provided with 14 coefficients for 53 elements . We need to extract the required data to calculate the specific heat , entropy and enthalpy for the required elements. Each elements will have a respective range of local temperatures.

The NASA polynomials have the form:

    Cp/R = a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4 -> Specific heat
    
    H/RT = a1 + a2 T /2 + a3 T^2 /3 + a4 T^3 /4 + a5 T^4 /5 + a6/T -> Enthalpy
    
    S/R  = a1 lnT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7 ->Entropy 

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 quantity of heat required to raise the temperature of one gram of a substance by one Celsius degree. The units of specific heat are usually calories or joules per gram per Celsius degree.

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

Function for specific heat :

% Specific heat calculation functionmm
function Cp = Specificheat (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) + (a18* (Temp).^2) + (a11* (Temp).^3) + (a12* (Temp).^4));
    
else 
    Cp = R*(al+ (a2*Temp) + (a3* (Temp).^2) + (a4* (Temp).^3) + (a5*(Temp).^14));
    % Else maximum coefficients are taken
end

Entropy :

The entropy of an object is a measure of the amount of energy which is unavailable to do work. Entropy is also a measure of the number of possible arrangements the atoms in a system can have. In this sense, entropy is a measure of uncertainty or randomness.

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

Function for 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
% Else maximum coefficients are taken
S = R.*(a1.*log(Temp) + (a2.*Temp)) + (a3*((Temp).^2)/2) + (a4*((Temp).^13)/3) + (a5*((Temp).^4)/4) + a7;
end

Enthalpy:

the sum of the internal energy and the product of the pressure and volume of a thermodynamic system.

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

Function for Enthalpy:

%Enthalpy Calculation function
function H = Enthalpy(al,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 maximum coefficients are taken
H = R*Temp.* (a1+((a2*Temp)/2) + ((a3*(Temp).^2)/3) + ((a4* (Temp).^3)/4) + ((a5* (Temp).^4)/5) +(a6./Temp));

end

Molecular weight calculation Function : There is a for loop iterated to find out the respective molecular weights for the species , this function is called at the final results collation process of the main program.

%function for calculation of molecular weight 

function Molecularweight = Molecular_weight(Species_SP)
% Standard Atomic Name
atomic_name = ['H', 'c','o', 'N', 'A'];
% Hydrogen, Carbon, Oxygen, Nitrogen, Argon
fprintf('Molecular weight of ');
fprintf(Species_SP);fprintf(' : ');

% Standard Atomic Number

Atomic Weight [1.008 12.011 15.999 14.007 39.948);

Starting Range of Molecule
Molecularweight = 0;
% Adding Atomic number to obtain particular value
for i= 1: length(Species_SP)
for j = 1:length(atomic_name)
if strcmpi(Species_SP(i), atomic_name(j))
Molecularweight = Molecularweight +Atomic_weight();
position j;
end
end

n = str2num(Species_SP(i));
if n>1
Molecularweight= Molecularweight + Atomic_weight (position)*(n-1);
end
end
fprintf('molecular weight %s', Species_SP);
fprintf('%d',Molecularweight);
disp(' '); 
end

Main Program :

clear all
close all
clc
% Make directory
mkdir 'C:\Users\vigne\Documents\MATLAB\File parsing'
% Reading Data file NASA Thermodynamic
% Open a file & 'r' Open the file for reading
f1 =fopen( 'THERMO.dat','r');
gt=fgetl(f1);

% Obtaing Temperatures - by Reading Header Information
global_temp=fgetl(f1);

%Spliting
t=strsplit(global_temp);

% Converting string cells to numbers using str2num,
global_low_temp=str2double(t{2});
global_medium_temp = str2double(t{3});
global_high_temp=str2double(t{4});

% Temperature range & Gas Constant
% Temperature Range assumed
Temp=linspace(global_low_temp, global_high_temp, 1000);

% Universal Gas Constant J/(mol-K)
R = 8.314;%

% skipping three comment lines in the file
for i=1:3
gt=fgetl(f1);
end
for i = 1:53
gt=fgetl(f1);
A=strsplit(gt,' ');
Species_SP=(A{1});
line_1=fgetl(f1);
a=strfind(line_1, 'E');

% Extracting higher and lower temperature coefficients
% Extracting FIRST 7 coefficients i.e HIGH- temperature coefficients

a1=line_1(1:a(1)+3);
a1=str2double(a1);           %1

a2=line_1(a(1)+4:a(2)+3);
a2=str2double(a2);          %2

a3=line_1(a(2) +4:a(3)+3);
a3=str2double(a3);            %3
a4=line_1(a(3) +4:a(4)+3);
a4=str2double(a4);           %4

a5=line_1(a(4)+4:a(5)+3);
a5=str2double(a5);           %5

line_2=fgetl(f1);
b=strfind(line_1, 'E');

a6=line_2(1:b(1)+3);
a6=str2double(a6);       %6

a7=line_2(a(1)+4:b(2)+3);
a7=str2double(a7);   %7


% Extracting SECOND Z coefficients - LOW-temperature coefficients
a8=line_2(b(2)+4:b(3)+3);
a8=str2double(a8);       %8

a9=line_2(b(3) +4:b(4)+3);
a9=str2double(a9) ;              %9

a10=line_2 (b(4)+4:b(5)+3);
a10=str2double(a10) ;             %10

line_3=fgetl(f1);
C=strfind(line_1, 'E'); 

a11=line_3(1:C(1)+3);
a11=str2double(a11);          %a11

a12=line_3(C(1)+4:C(2)+3);     %12
a12=str2double( a12);

a13=line_3(C(2)+4:C(3)+3);
a13=str2double(a13);  %13

a14=line_3(C(3)+4:C(4)+3);
a14=str2double(a14);   %14

% Specific heat calculation
Cp = Specificheat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp);

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

% Entropy calculation

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

%Molecular weight calculation 

Molecularweight = Molecular_weight(Species_SP);

% Read & Write Molecular list
Molecular_list = fopen('Molecular weight.txt', 'a+');
fprintf(Molecular_list,' Molecular weight of %s is %d n', Species_SP, Molecularweight);
fclose(Molecular_list);

% plotting(CP,H,S, Temp, Species_SP)

% Saving as a folder
cur_dir=pwd;
mkdir(Species_SP);
cd(Species_SP);

% plot 1, Specific heat (CP) vs. Temperature Range (Temp)
figure(1)
plot(Temp,Cp, 'linewidth', 2,'color', 'r');
xlabel('Temperature [K]');
ylabel ('Specific Heat [k]');
title( sprintf("varying specific Heat with Temperature Range of %s",Species_SP));
grid on
saveas(figure(1),'Specific Heat.jpg');

% Plot 2, Enthalpy (H) vs Temperature Range (Temp)
figure(2)
plot(Temp,H,'linewidth',2,'color','G');
xlabel('Temperature [K]');
ylabel('Enthalpy in [kJ/(mol-K)]');

title(fprintf('Varying Enthalpy with Temperature range of %s' ,Species_SP));

grid on
saveas(figure(2), 'Ethalpy.jpg');

%Plot 3, Entropy (s) vs Temperature Range (Temp)
figure(3)
plot(Temp,S, 'linewidth',2, 'color', 'b');
xlabel('Temperature [K]');
ylabel('Entropy in [k]/(mol-k)]');
title(sprintf( 'Varying Entropy with Temperature Range of %s', Species_SP));
grid on
saveas(figure(3), 'Entropy.jpg');
cd(cur_dir);
% Molecular weight
M = Molecular_weight(Species_SP);
f2= fopen('molecular mass', 'w');
fprintf('%s n % n %f n', 'species', 'Molecular mass',M)
fclose(f2);
end

Explanation :

• The program clearly contains the command lines with the execution flow . Please refer the commands for clear understanding
• Each required physical quantity is found thro by calling the respective functions .
• Firstly, we are trying to plot the graph for Specific heat (CP) vs. Temperature Range (Temp) . The graph is obtained for each expected species .
• Similar way we are trying to plot graphs for Enthalpy (H) vs Temperature Range (Temp) and Entropy (s) vs Temperature Range (Temp) .
• We are finally trying to calculate the molecular weight for the species by calling the function declared previopsuly for molecular weight.

Output :

1.Directory Snapshot

2.Molecular weight [Command Window]

Sample output plots :

Species: O2

Species : CO2

Errors :

• Got many suntax errors [screenshot was not taken]
• Functions declaration was having issues .
• The specific heat graph obtained was -ve , checked the function declared for the specific heat there was a '-' sign implemented by mistake.

References:

1. https://www.youtube.com/watch?v=xH2ArJK2tO8
2. https://in.mathworks.com/matlabcentral/answers/497882-parsing-text-files-for-beginners
3. https://in.mathworks.com/help/matlab/function-basics.html

***********************************Thank You********************************************

Vigneshwar Viswanathan