TO PARSE THE ‘NASA THERMODYNAMIC DATA’ USING MATLAB/OCTAVE
AIM: TO PARSE THE ‘NASA THERMODYNAMIC DATA’ AND TO CALCULATE DIFFERENT THERMODYNAMIC PROPERTIES USING MATLAB/OCTAVE.
OBJECTIVE:
- To read the Nasa thermodynamic data file.
- To write a function that extracts the 14 coefficients and calculates the specific heat, enthalpy, entropy for all species in the data file.
- To calculate the molecular weight of each species and to display it in the command window.
- To plot the graphs of specific heat, enthalpy, entropy with temperature range.
- To save the plots as images with appropriate names and to dump them in separate folders for each species with suitable name using program.
OVERVIEW:
Parsing in computer languages refers to syntactic analysis of the input code into its components parts in order to facilitate the writing of compilers and interpreters. It can be used to describe split or separation. Parsing a file means reading in a data stream of some sort and building an in-memory model of the semantic content of that data. This aims at facilitating performing some kind of transformation on the data.
Where,
a1, a2, a3, a4, a5, a6, and a7 are the numerical coefficients provided in NASA thermodynamic files.
R is the universal gas constant = 8.314 J/mol. K
T is the local temperature
BODY OF CONTENT:
1. Program for Specific Heat -
%TO CALCULATE SPECIFIC HEAT
function c_p = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,T,global_min_temp,global_mid_temp,global_max_temp)
if T>global_mid_temp
c_p = R*(a1 + a2.*T + a3.*(T.^2) + a4.*(T.^3) + a5.*(T.^4));
else
c_p = R*(a8 + a9.*T + a10.*(T.^2) + a11.*(T.^3) + a12.*(T.^4));
endEXPLANATION OF CODE:
- This function is used to calculate the Specific heat for each species.
- a1, a2, a3, a4, a5, a6, a7 are numerical coefficients for ‘high temperature range’.
- a8, a9, a10, a11, a12, a13, a14 are numerical coefficients for ‘Low temperature range’.
- ‘If-else’ conditional statement is used to assign condition of different temperature range.
2. Program for Enthalpy -
%TO CALCULATE ENTHALPY
function H = enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,T,global_min_temp,global_mid_temp,global_max_temp)
if T>global_mid_temp
H = R.*T.*(a1 + a2.*(T./2) + a3.*((T.^2)./3) + a4.*((T.^3)./4) + a5.*((T.^4)./5) + a6./T);
else
H = R.*T.*(a8 + a9.*(T./2) + a10.*((T.^2)./3) + a11.*((T.^3)./4) + a12.*((T.^4)./5) + a13./T) ;
endEXPLANATION OF CODE:
- This function is used to calculate the Enthalpy for each species.
- a1, a2, a3, a4, a5, a6, a7 are numerical coefficients for ‘high temperature range’.
- a8, a9, a10, a11, a12, a13, a14 are numerical coefficients for ‘Low temperature range’.
- ‘If-else’ conditional statement is used to assign condition of different temperature range.
3. Program for Entropy -
%TO CALCULATE ENTROPY
function S = entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,T,global_min_temp,global_mid_temp,global_max_temp)
if T>global_mid_temp
S = R*(a1.*log(T) + a2.*T + a3.*(T.^2)./2 + a4.*(T.^3)./3 + a5.*(T.^4)./4 + a7);
else
S= R*(a8.*log(T) + a9.*T + a10.*(T.^2)./2 + a11.*(T.^3)./3 + a12.*(T.^4)./5 + a14);
end EXPLANATION OF CODE:
- This function is used to calculate the Entropy for each species.
- a1, a2, a3, a4, a5, a6, a7 are numerical coefficients for ‘high temperature range’.
- a8, a9, a10, a11, a12, a13, a14 are numerical coefficients for ‘Low temperature range’.
- ‘If-else’ conditional statement is used to assign condition of different temperature range.
4. Program for plots -
% Plot for all Species
function plotting_properties(c_p,H,S,T,element_name)
% Creating folder
cd(pwd);
mkdir(element_name);
cd(element_name);
% Plot-1 "SPECIFIC HEAT vs TEMPERATURE"
figure(1)
plot(T,c_p,'linewidth',3,'color','g');
xlabel('Temperature [K]');
ylabel('Specific Heat [KJ]');
title(sprintf('Specific Heat vs Temperature range for %s:',element_name));
grid on
saveas(figure(1),'Specific Heat vs Temperature.jpg');
% Plot-2 "ENTHALPY vs Temperature"
figure(2)
plot(T,H,'linewidth',3,'color','b');
xlabel('Temperature [K]');
ylabel('Enthalpy [KJ]');
title(sprintf('Enthalpy vs Temperature range %s:',element_name));
grid on
saveas(figure(2),'Enthalpy vs Temperature.jpg');
% Plot-3 "ENTROPY vs TEMPERATURE"
figure(3)
plot(T,S,'linewidth',3,'color','r');
xlabel('Temperature[K]');
ylabel('Entropy [KJ]');
title(sprintf('Enthalpy vs Temperature range %s:',element_name));
saveas(figure(3),'Entropy vs Temperature.jpg');
cd .. ;
endfunction
EXPLANATION OF CODE:
- This function is used to plot the graphs for each species.
- ‘cd’ command is used to change the working dictionary.
- ‘pwd’ command is to identify the current dictionary.
- ‘mkdir’ command is used to make new dictionary.
- ‘title’ command is used to entitle the graph with appropriate species name.
- ‘saveas’ command is used to save the plots as image file.
5. Program to calculate molecular weight -
% TO CALCULATE MOLECULAR WEIGHT
function [MW] = Molecular_weight_cal(element_name)
elements = ['H' 'C' 'N' 'O' 'Ar'];
% ATOMIC WEIGHT FOR EACH ELEMENT
atomic_wt = [1 12 14 16 40];
MW = 0;
for i = 1:length(element_name)
for j = 1:length(elements)
if strcmp(element_name(i), elements(j))
MW = MW + atomic_wt(j);
% STORING THE POSITION OF MATCHED ELEMENT
p = j;
end
end
n = str2num(element_name(i));
if n > 1
MW = MW + atomic_wt(p)*(n-1);
end
end
endEXPLANATION OF CODE:
- This function is used to calculate the molecular weight for each species.
- Atomic weight is assigned to the elements. [ H > 1, C>12, N>14, O>16, Ar>40]
- ‘strcmp’ command is used to compare the name of the element with the string array.
- ‘str2num’ command is used to convert character array or string to numeric array.
6. Program to write molecular weight to a file -
% TO WRITe MOLECULAR WEIGHT OF ELEMENTS TO A FILE
function write_file(element_name,Molecularweight)
mol_list = fopen('mol_wt.dat','a+');
fprintf(mol_list, 'Molecular wieght of %s \n',element_name,Molecularweight);
fclose(mol_list);
endEXPLANATION OF CODE:
- A new file is created with name 'mol_wt.dat'and it is opened in append mode denoted by 'a+' so that we can write the molecular weight of all elements in a single file.
- ‘fopen’ is used to open the file named ‘mol_wt.dat’.
- Then ‘fprintf’ is used to write the molecular weight data with element name is the ‘mol_wt.dat’ file.
Final Program -
CODE:
clear all
close all
clc
R = 8.314; % Universal gas constant in [J/mol.K]
% Read the data file
f1 = fopen('THERMO.dat','r'); % r is for read mode
% For first line
L1 = fgetl(f1);
% For second line
% To get the global min, med, max temperatures
L2 = fgetl(f1);
A = strsplit(L2);
global_min_temp = str2num(A{2});
global_mid_temp = str2num(A{3});
global_max_temp = str2num(A{4});
% Temperature range
T = linspace(global_min_temp,global_max_temp,500);
% To skip the unwanted lines
for i = 1:3
fgetl(f1);
end
for i = 1:53
I = fgetl(f1); % Information
element_name = strtok(I);
% line-1,2 are required for all high temp. coeffs.
line_1 = fgetl(f1);
line_2 = fgetl(f1);
% To find the position of "E"
a = strfind(line_1,'E');
b = strfind(line_2,'E');
% High temperature Coefficients
a1 = str2num(line_1(1:a(1)+3));
a2 = str2num(line_1(a(1)+3:a(2)+3));
a3 = str2num(line_1(a(2)+4:a(3)+3));
a4 = str2num(line_1(a(3)+4:a(4)+3));
a5 = str2num(line_1(a(4)+4:a(5)+3));
a6 = str2num(line_1(1:b(1)+3));
a7 = str2num(line_2(b(1)+4:b(2)+3));
% line-3 is also required for low temp. coeffs.
line_3 = fgetl(f1);
c = strfind(line_3,'E');
% Low temperature Coefficients
a8 = str2num(line_2(b(2)+4:b(3)+3));
a9 = str2num(line_2(b(3)+4:b(4)+3));
a10 = str2num(line_2(b(4)+4:b(5)+3));
a11 = str2num(line_3(1:c(1)+3));
a12 = str2num(line_3(c(1)+4:c(2)+3));
a13 = str2num(line_3(c(2)+4:c(3)+3));
a14 = str2num(line_3(c(3)+4:c(4)+3));
% To calculate Specific heat
c_p = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,T,global_min_temp,global_mid_temp,global_max_temp);
% To calculate Enthalpy
H = enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,T,global_min_temp,global_mid_temp,global_max_temp);
% To calculate Entropy
S = entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,T,global_min_temp,global_mid_temp,global_max_temp);
% Plots
plotting_properties(c_p,H,S,T,element_name);
% molecular weight
mol_wt(i) = Molecular_weight_cal(element_name);
disp(mol_wt(i));
Molecular_weight = strcat(mol_wt(i),element_name)
% write to file [MW]
Molecular_weight = mol_wt(i);
write_file(element_name,Molecular_weight);
end
fclose(f1)EXPLANATION OF CODE:
- Clear all, close all, clc command is used to clear all previous commands and plots.
- ‘fopen’ is used to open the file named ‘THERMO.dat’.
- ‘r’ is used for read mode.
- ‘fgetl’ is to read the file line-by-line.
- ‘strsplit’ command is for splitting the string or character vector at specified delimiter.
- ‘str2num’ command is used to convert character array or string to numeric array.
- ‘linspace’ is used to generate the temperature range.
- To read all the 53 species one by one and plot them is different folder ‘for-loop’ is used.
- ‘strfind is to find one string within another.
- All functions defined previously are written one by one.
- ‘fclose’ is used to close the data file.
ERROR:
1. Error in enthalpy code, [forgot one ‘)’] –
2. Error in entropy code, [assigned scaler dot product] –
3. Error as ‘molecular_weight’ undefined [cd command was absent at the end of the ‘plotting_properties’ function] -
RESULTS:
SCREENSHOT FOR DICTIONARY FOLDER -
SCREENSHOTS OF MOLECULAR WEIGHT -
ALL 14 COEFFICIENTS [a1, a2, a3, a4, a5, a6, a7 are for 'High temperature range' and a8, a9, a10, a11, a12, a13, a14 are for 'Low temperature range'] -
PLOTS
1. Plots for O2 –
2. Plots for N2 –
3. Plots for CO2 –
REFERENCES: