File Parsing of NASA s Thermodynamic data Using Matlab
AIM:- To Parse the NASA Thermodynamic Properties Data file and introduce the Required properties.
OBJECTIVES:-
- To Extract the 14 coefficients and to calculate the Enthalpy, Entropy, and Specific heats for all the species in the data file.
- Calculate the molecular weight of each species.
- Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species.
INTRODUCTION:-
NASA\'s Thermodynamic properties file consists of multiple species of gases, their respective operating temperature ranges, and coefficients to calculate thermodynamic properties. There are 14 coefficients for each specie of which the first 7 are \"High Temperature Coefficients\" and the next 7 are \"low temperature coefficients\". 14 Coefficients are split into 3 lines, 5 each in the first 2 and 4 in the 3rd. This makes reading the line and isolating the coefficients a challenge, so understanding the positions of coefficients and checking for repitition across both lines is the key here.
Once Extracted, the below equations are used to calculate Specific heat(Cp), Enthalpy(H),and Entropy(s) for Each of the 53 species in the file.
The values generated are then plotted vs Local temperatures.
EXPLANATION:-
- Code starts with the clear all, close all and clc. Then value of universal gas constant \'R\' is provided.
- In order to open the file, \'fopen\' command is used, \'THERMO.dat\' is the file namew \'r\' relates to read a file.
- \'fgetl\' command is used to get the first two lines l1 &l2. l1 gives the name of the file and l2 gives the temperature points which is global minimum, global medium & global maximum temperature. \'strsplit\' command is used to split the string by default set of sequence and gather the corresponding temperature data points.
- \'str2num\' was initially used to convert the string to number. However str2double is used to convert the strings into precise value.
- Linspace command is used to get the 400 values of the temperature that varies between global minimum temperature and global maximum temperature.
- After the temperatures, the next 3 lines from the file are not required.
- Then a for loop is used for the 53 constituents mentioned in the data file to gather the name of the compound and is given to the species_name from the fgetl and strsplit commands.
- d_line1,dline2,& d_line3, these 3 lines contains the 14 coefficients for determining the properties.
- strfind command is used to identify the letter E from the strings.
- htc1,htc2,htc3,htc4,htc5,htc6,htc6,htc7 represents the high temperature coefficients.
- ltc1,ltc2,ltc3,ltc4,ltc5,ltc6,ltc7 represents the low temperature coefficients.
- The strings are converted to values using the str2double and by identifying E in the strings accordingly.
- mkdir command is used to make the directory for all the values and plots are determined.
- Then Separate function is used to calculate Enthalpy, Entropy & Specific heat.
- Here global mid temperature, high, & low coefficients,Universal gas constant re used as an input.
- Then, if the temperature is equal or below the global medium temperature, the Enthalpy should be calculated by using low temperature coefficients,Otherwise high temperature coefficients are used to calculate Enthalpy. This is achieved by using \' if-else \' command.
- Then, if the temperature is equal or below the global medium temperature, the Enthalpy should be calculated by using low temperature coefficients,Otherwise high temperature coefficients are used to calculate Entropy. This is achieved by using \' if-else \' command.
- Then, if the temperature is equal or below the global medium temperature, the Enthalpy should be calculated by using low temperature coefficients,Otherwise high temperature coefficients are used to calculate Specific heat. This is achieved by using \' if-else \' command.
MAIN CODE:-
clear all
close all
clc
R = 8.314; %Universal Gas Constant
f1 = fopen(\'THERMO.dat\', \'r\')
a=fgetl(f1)
%Reading first line
fprintf(\'Data file title: %s\',a)
%Reading Global Temperatures line
% Splitting and storing into cell array, & Converting to number
b=fgetl(f1)
globaltemps = strsplit(b,\' \')
globallow = str2double(globaltemps{2})
globalmid = str2double(globaltemps{3})
globalhigh = str2double(globaltemps{4})
%Counting total number of lines traversed so far
count = 2;
%Traversing commented lines
for i = 1:3
fgetl(f1);
count = count +1;
end
%Lines is the total number of lines in the document that we will use in the loop later on
lines = 0;
f2 = fopen(\'THERMO.dat\', \'r\');
while fgetl(f2) ~= -1
lines = lines +1;
end
%Below is the total number of lines after the commented lines in the DAT file
lines = lines - count;
fclose(f2);
% Variable Extraction and Calculation
for i = 1:lines/4 %Every 4th line is an element\'s start, so we skip 4 lines in the counting
% Starting off with the species line
sline = strsplit(fgetl(f1),\' \');
fprintf(\'Calculation for: %s\', sline{1})
%Extracting local Low, Medium, High Temperatures
locallow = str2double(sline{length(sline)-3})
localmed = str2double(sline{length(sline)-1})
localhigh = str2double(sline{length(sline)-2})
% Local temperature range
T = linspace(locallow, localhigh, 200);
% Molecular weight calculation
% Mol.Wt. Calaculated Separately in function \"molwtcalc2\"
fprintf(\'Molecular weight of %s is %d\', sline{1}, molwtcalc2(sline(1)))
% Coefficient Extraction starts here
dline1 = fgetl(f1);
a = strfind(dline1, \'E\') % Changed to strfind from findstr later as it is faster
dline2 = fgetl(f1);
%b = findstr(dline2, \'E\') % Since position of coefficients is the same in all lines, one line is enough
dline3 = fgetl(f1);
%c = findstr(dline3, \'E\')
% High Temperature Coefficients
htc1 =str2double(dline1(1:a(1)+3));
htc2 =str2double(dline1(a(1)+4:a(2)+3));
htc3 =str2double(dline1(a(2)+4:a(3)+3));
htc4 = str2double(dline1(a(3)+4:a(4)+3));
htc5 = str2double(dline1(a(4)+4:a(5)+3));
htc6 = str2double(dline2(1:a(1)+3));
htc7 = str2double(dline2(a(1)+4:a(2)+3));
% Low Temperature Coefficients
ltc1 = str2double(dline2(a(2)+4:a(3)+3));
ltc2 = str2double(dline2(a(3)+4:a(4)+3));
ltc3 = str2double(dline2(a(4)+4:a(5)+3));
ltc4 = str2double(dline3(1:a(1)+3));
ltc5 = str2double(dline3(a(1)+4:a(2)+3));
ltc6 = str2double(dline3(a(2)+4:a(3)+3));
ltc7 = str2double(dline3(a(3)+4:a(4)+3));
% Calculation of Cp, H and S
for k = 1:length(T)
if T(k) > localmed %high temp coeff used for temps > median temperature
Cp(k) = R*(htc1 + (htc2*T(k)) + (htc3*(T(k)^2)) + (htc4*(T(k)^3)) + (htc5*(T(k)^4)));
H(k) = R*(htc1 + ((htc2*T(k))/2) + ((htc3*(T(k))^2)/3) + ((htc4*(T(k))^3)/4) + ((htc5*(T(k))^4)/5) + (htc6/T(k)));
S(k) = R*(htc1*log(T(k))+ (htc2*T(k)) + ((htc3*(T(k)^2))/2) + ((htc4*(T(k)^3))/3) + ((htc5*(T(k)^4))/4) + htc7);
else % Low temp coefficients used for Temps < median temperature
Cp(k) = R*(ltc1 + (ltc2*T(k))+ (ltc3*(T(k)^2)) + (ltc4*(T(k)^3)) + (ltc5*(T(k)^4)));
H(k) = R*(ltc1 + ((ltc2*T(k))/2) + ((ltc3*(T(k))^2)/3) + ((ltc4*(T(k))^3)/4) + ((ltc5*(T(k))^4)/5) + (ltc6/T(k)));
S(k) = R*(ltc1*log(T(k)) + (ltc2*T(k)) + ((ltc3*(T(k)^2))/2) + ((ltc4*(T(k)^3))/3) + ((ltc5*(T(k)^4))/5) + ltc7);
end
end
% Creating directories
%mkdir([\'Parsing_Nasa_thermodynamic_data\',sline{1}]) %Stored in separate directory
% cd([\'Parsing_Nasa_thermodynamic_data\',sline{1}])
% Plotting specific Heat, (Cp vs T)
% a = i; Earlier \'i\' is used here to place in \"figure(i)\" to dispaly
% different no. of figures
figure(1)
plot(T, Cp, \'linewidth\', 3, \'color\', \'r\')
xlabel(\'Temperature(k)\')
ylabel(\'Specific Heat at Constant Pressure (J/KgK)\')
title([\'Specific Heat vs Temperature for %s\', sline{1}])
saveas(1, ([\'Specific Heat vs Temperature for %s\', sline{1}]), \'jpeg\')
% Plotting Enthalpy, (H vs T)
% b = i+53;
figure(2)
plot(T, H, \'linewidth\', 3, \'color\', \'g\')
xlabel(\'Temperature(k)\')
ylabel(\'Enthalpy (J)\')
title([\'Enthalpy vs Temperature for %s \', sline{1}])
saveas(2, ([\'Enthalpy vs Temperature for %s \', sline{1}]), \'jpeg\')
% Plotting Entropy, (S vs T)
% c = b+53;
figure(3)
plot(T, S, \'linewidth\', 3, \'color\', \'b\')
xlabel(\'Temperature(k)\')
ylabel(\'Entropy (J/K)\')
title([\'Entropy vs Temperature for %s \', sline{1}])
saveas(3, ([\'Entropy vs Temperature for %s\', sline{1}]), \'jpeg\')
cd ..
disp(\'_______________________\')
cd(\'C:\\c data\\Desktop\\Matlab codes\')
% Output window separator line for better viewing
end
fclose(f1);
CODE FOR MOLECULAR WEIGHT:-
function [molwt] = molwtcal2(specie)
molwt = 0;
elements = split(specie,\"\");
for j = 2:length(elements)-1
n = str2num(elements{j+1});
if elements{j} == \'C\'
if n >=1
molwt = molwt + (n*12);
else
molwt = molwt + 12;
end
elseif elements{j} == \'H\'
if n >= 1
molwt = molwt + (n*1);
else
molwt = molwt + 1;
end
elseif elements{j} == \'O\'
if n >=1
molwt = molwt + (n*16);
else
molwt = molwt + 16;
end
elseif elements{j} == \'N\'
if n >=1
molwt = molwt + (n*14);
else
molwt = molwt + 14;
end
elseif (elements{j} == \'A\') && (elements{j+1} == \'R\') %AR = \'A\' and \'R\' in 2 separate cell locations
if n >=1
molwt = molwt + (n*40);
else
molwt = molwt +40;
end
end
end
molwt;
end
Plots:-
1) O2:-
2) CO2:-
3) N2:-
WINDOW SCREENSHOTS:-
WORKSPACE:-
OUTPUT WINDOW:-
REFERENCES:-
- https://in.mathworks.com/help/gads/ga.html#mw_f69991d0-5b23-4023-8c8e-e2890f2474d5
- https://en.wikipedia.org/wiki/MathWorks
- https://drive.google.com/file/d/1AikCHDlz9s_fu8whE81qee_GD4AyYpCP/view?usp=sharing.