Parsing NASA’s Thermodynamic data and evaluating properties such as Cp, H, S.
Aim:
To parse NASA’s Thermodynamic data and evaluate properties such as Cp, H, S.
Objective:
- To parse the thermodynamic data and save it.
- To calculate the thermodynamic properties for each element.
- To plot the thermodynamic properties.
- To calculate the molecular weight of each element.
Governing Equations:
NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S. The polynomials are shown below:
Where R is the universal gas constant = 8.314 J/mole-K.
The first 7 coefficients are high temperature coefficients and the second 7 coefficients are low temperature coefficients.
Workflow:
Program Code:
clear all
close all
clc
% saving all lines of text in a structure array
ct = 1;
thermo = fopen('THERMO.dat','r'); % opening the data file
tline = fgetl(thermo); % skipping first line
while ischar(tline)
tline = fgetl(thermo); % loop to extract all lines one by one
x(ct).data = tline; % line is stored in a structure array
ct = ct+1;
end
% getting global temperatures
gt1 = x(1).data;
gt2 = split(gt1); % splitting the string
for i = 2:length(gt2)
global_temp(i-1) = str2double(gt2{i}); % converting temperature strings to numbers
end
% getting local temperatures for all elements
k = 1;
for j=5:4:ct
lt1 = x(j).data;
lt2 = split(lt1);
for i = numel(lt2)-3 : numel(lt2)-1
if numel(lt2) == 12
T(k,i-8) = str2double(lt2{i}); % converting temperature strings to numbers
elseif numel(lt2) == 11 % using different if else satements based
T(k,i-7) = str2double(lt2{i}); % on different number of strings in the line
elseif numel(lt2) == 10
T(k,i-6) = str2double(lt2{i});
elseif numel(lt2) == 9
T(k,i-5) = str2double(lt2{i});
elseif numel(lt2) == 8
T(k,i-4) = str2double(lt2{i});
end
end
k = k+1;
end
% getting coefficients for all elements
k = 1;
for i = 6:4:ct-6
% getting coefficients from 1st line
tline = x(i).data;
b = strfind(tline,'E'); % finding the character E and using it to split the string
a(k,1) = str2double(tline(1:b(1)+3));
a(k,2) = str2double(tline(b(1)+4:b(2)+3));
a(k,3) = str2double(tline(b(2)+4:b(3)+3));
a(k,4) = str2double(tline(b(3)+4:b(4)+3));
a(k,5) = str2double(tline(b(4)+4:b(5)+3));
% geting coefficients from 2nd line
tline = x(i+1).data;
b = strfind(tline,'E');
a(k,6) = str2double(tline(1:b(1)+3));
a(k,7) = str2double(tline(b(1)+4:b(2)+3));
a(k,8) = str2double(tline(b(2)+4:b(3)+3));
a(k,9) = str2double(tline(b(3)+4:b(4)+3));
a(k,10) = str2double(tline(b(4)+4:b(5)+3));
% geting coefficients from 3rd line
tline = x(i+2).data;
b = strfind(tline,'E');
a(k,11) = str2double(tline(1:b(1)+3));
a(k,12) = str2double(tline(b(1)+4:b(2)+3));
a(k,13) = str2double(tline(b(2)+4:b(3)+3));
a(k,14) = str2double(tline(b(3)+4:b(4)+3));
k = k+1;
end
% for getting names of all elements
k = 1;
elements = cell(1,length(a)); % making a cell array to store the element names
for i = 5:4:ct-7
elements{k} = strtok(x(i).data); % saving the first string using stringtok
k = k+1;
end
% calculating Cp, H and S
for j = 1:length(a)
k = 1;
for i = T(j,1):T(j,3) % for lower temperature range
cp_low(j,k) = 8.314*(a(j,8) + a(j,9)*i + a(j,10)*i^2 + a(j,11)*i^3 + a(j,12)*i^4);
h_low(j,k) = 8.314*i*(a(j,8) + (a(j,9)*i)/2 + (a(j,10)*i^2)/3 + (a(j,11)*i^3)/4 + (a(j,12)*i^4)/5 + a(j,13)/i);
s_low(j,k) = 8.314*(a(j,8)*log(i) + a(j,9)*i + (a(j,10)*i^2)/2 + (a(j,11)*i^3)/3 + (a(j,12)*i^4)/4 + a(j,14));
k = k+1;
end
k = 1;
for i = T(j,3):T(j,2) % for higher temperature range
cp_high(j,k) = 8.14*(a(j,1) + a(j,2)*i + a(j,3)*i^2 + a(j,4)*i^3 + a(j,5)*i^4);
h_high(j,k) = 8.14*i*(a(j,1) + (a(j,2)*i)/2 + (a(j,3)*i^2)/3 + (a(j,4)*i^3)/4 + (a(j,5)*i^4)/5 + a(j,6)/i);
s_high(j,k) = 8.14*(a(j,1)*log(i) + a(j,2)*i + (a(j,3)*i^2)/2 + (a(j,4)*i^3)/3 + (a(j,5)*i^4)/4 + a(j,7));
k = k+1;
end
end
% calculating molecular weights
atoms = ['H','C','O','N','A','S']; % making sring of elements with corresponding atomic weight
atomic_wt = [1 12 16 14 40 32];
for i = 1:length(elements)
mw(i) = 0; % initialising mw to 0
element = elements{i};
for j = 1:length(element)
for k = 1:length(atoms)
if strcmp(element(j),atoms(k)) % comparing string to known element
mw(i) = mw(i) + atomic_wt(k);
y = k;
end
end
n = str2double(element(j));
if n>1
mw(i) = mw(i) + atomic_wt(y)*(n-1); % multiplying if element has multiple of same atom
end
end
end
% displaying molecular weights
for i = 1:length(a)
fprintf('The molecular weight of %s is = %d \n',elements{i},mw(i));
end
% plotting graphs and saving in folders
mkdir('NASA_Graphs'); % making new parent directory
cd('NASA_Graphs'); % changing current directory
for i = 1:length(a)
t_low = T(i,3) - T(i,1) + 1; % temporarily storing temperature range
t_high = T(i,2) - T(i,3) + 1;
mkdir(elements{i}); % making directory for each element
cd(elements{i});
figure(1);
hold on
plot(T(i,1):T(i,3),cp_low(i,1:t_low),'linewidth',2); % plotting the curve
plot(T(i,3):T(i,2),cp_high(i,1:t_high),'linewidth',2);
title('Cp vs T');
xlabel('Temperature');
ylabel('Specific Heat');
legend('Low Temperature Range','High Temperature Range');
hold off
saveas(gcf,'Cp_vs_T.png'); % saving figure
close % closing figure to prevent overlap
figure(2);
hold on
plot(T(i,1):T(i,3),h_low(i,1:t_low),'linewidth',2);
plot(T(i,3):T(i,2),h_high(i,1:t_high),'linewidth',2);
title('H vs T');
xlabel('Temperature');
ylabel('Enthalpy');
legend('Low Temperature Range','High Temperature Range');
hold off
saveas(gcf,'H_vs_T.png');
close
figure(3);
hold on
plot(T(i,1):T(i,3),s_low(i,1:t_low),'linewidth',2);
plot(T(i,3):T(i,2),s_high(i,1:t_high),'linewidth',2);
title('S vs T');
xlabel('Temperature');
ylabel('Entropy');
legend('Low Temperature Range','High Temperature Range');
hold off
saveas(gcf,'S_vs_T.png');
close
cd .. % going back to parent directory for next iteration
endCode Explanation:
1. Clearing workspace and saving all lines of text from data file into a structure array using while loop.
2. Getting global temperatures from the structure array using while loop. Splitting the entire line into multiple strings and saving each temperature string after converting to a number using str2double.
3. Getting local temperatures from the structure array using concatenated for loops.
4. Recording every 4th line which contains then temperature values and using if elseif statements based on different possibilities of string sizes to prevent null entries.
5. Reading the polynomial coefficients from every 4th line using for loop.
6. Finding the alphabet ‘E’ in every line and saving each string with appropriate offset on both sides after converting to number using str2double.
7. Doing the same for the 2nd and 3rd line using i + 1 and i + 2.
8. Reading the names of all elements and storing them in a cell array.
9. strtok function records the first string before encountering a whitespace. We use it to record the element name which is written in the starting of every 4th line.
10. Calculating specific heat, enthalpy and entropy for both low and high local temperature range for each element using for loop by using the appropriate coefficients.
11. Calculating the molecular weights for each element.
12. Defining the reference atoms and their corresponding atomic weights.
13. Comparing each atom of each element with the reference elements and increasing the molecular weight accordingly.
14. Using if else statement to check for multiatomic elements and adjusting their weights accordingly.
15. Displaying the molecular weights one by one using for loop.
16. Plotting graphs of Cp vs T for all elements for both low and high temperature range.
17. Making a directory for storing the graphs using mkdir command and changing to that directory.
18. Saving the figure as a png in the appropriate element folder.
19. Doing the same for H vs T and S vs T.
Errors:
Only Syntax errors were encountered which were resolved accordingly.
Output:
The molecular weights were displayed as shown:
- Plots for O2 are as shown:
- Plots for CO2 are as shown:
- Plots for N2 are as shown;