Thermodynamic parameter calculation Octave code from NASA Coefficient file
Introduction
Thermodynamics is basically a science dealing with interconversion of heat and work. However, during this process, numerous state parameters change, in turn affecting the order of system, direction of process, state properties, kinectics of reaction. Apart from Pressure, Volume and temperature, enthalpy, entropy plays a vital role in calculating quality and quantity of heat and work interaction. These are species dependent and temperature dependent.
NASA, in late 1990s, compiled the exhaustive list of coefficients which could be used to calculate the values of specific heat at constant pressure, enthalpy and entropy. These coefficients are species specific. However the generalised equation to calculate the Specific heat at constant pressure, enthalpy and entropy remains the same for all species. However, for the same species the coefficients used is divided to two categories depending on temperature a1 to a7 are used for low temperature and a8 to a14 for high temperature where a1 to a14 represents coefficients of the species.
where cp stands for Specific heat at constant pressure, R stands for universal gas constant, T for temperature, H for enthalpy and S for entropy.
Reading the data from the dat file having over 300 species and calculating each parameter individually can take a lot of time and considerable human effort. Hence, this code is developed to overcome such monotonous repetitive tasks.
The code could be used in three different fashion with minor changes to the code
Mode 1: User requires specific species data and the plots for quick visualisation. The code scans the dat file and returns the values for the species user is interested in, plots the data for quick visualisation and saves the data to the directory for the future reference.
Mode 2: User doesnt input any specific species. But the code only takes in the Nasa dat file as the input and calculates values for each species in the NASA dat file and plots data for quick visualisation for each species and stores the plot for each species in the separate folders for future reference.
Mode 3: This is the hybrid of Mode 1 and Mode 2. User requires specific species data, but the code scans the dat file, calculates values for each species in NASA dat file and plots data for quick visualisation for each species and stores the plot in the separate folders and returns the plots the user is interested in from the folders in the directory.
In this present work, Mode 1 is enabled and presented. However, mode 2 was used for generating plots and data for all the species in the dat file. Transition from mode 1 to mode 2 can be easily done by commenting few lines as mentioned in the code.
Octave code:
The logic of code for mode 1 is as follows;
1) Code does the intial cleanup
2) Takes user input for species, user is interested in, inputs are made case insensitive
%This is coded by Sweekar, [email protected], please contact for queries
%This is an intellectual property of Sweekar, conatct before copying
%This code is used to read from NASA coefficients file for different species
%The code allows user to enter the desired species for which the calculation required
%User is expected to enter the species name by its chemical convention
%Code gets the coefficients from the NASA file and determines enthalpy, entropy and specific heat
%Plots are created for easy visualisation of calculation
%This code will take user inputs, generate plots, calculate molecular weight of species, create the directory and store the plots
%A minor change to this code will enable in generating plots, directories, molecular weight for all the species in Thermo.dat file
%cleanup
close all
clear all
clc
R = 8.314; %Universal gas constant
num_species = 53; %Number of species in the NASA dat file
%Input from user
fprintf("..............NASA Coefficients and Properties calc...............\n")
fprintf("\n The NASA file referenced has %d species\n",num_species)
fprintf("DISCLAIMER: If no plot is obtained, then species is not listed in dat file\n")
fprintf("Please change the details and specify the number of species accurately\n")
fprintf("Plesae enter the desired species name by its chemical formula\n") %Comment this line if program is used to generate properties of all species in dat file
fprintf("Enter within quotes, ex, if oxygen is required enter 'o2'\n")%Comment this line if program is used to generate properties of all species in dat file
User_species = input("Enter species of interest within single quotes: ");%Comment this line if program is used to generate properties of all species in dat file
User_species = upper(User_species); %Making case insensitive
nasa = fopen('THERMO.dat','r'); %Reading Nasa file
fgetl(nasa);3) Opens the NASA dat file and reads headings, stores Global temperatures.
4) A while loop is created with combination a Logical AND operator, Switch variable and Counter variable
5) Within the loop, if-else condition is utilised to compare the species user is interested in with the species header line.
Loctemp = fgetl(nasa);
B = strsplit(Loctemp,' ');
Loc_species = B{1};
w = strcmp(User_species,Loc_species); %Comment this line if program is used to generate properties of all species in dat file
if(w==1) %Comment this line if program is used to generate properties of all species in dat file
loop = 0;6) When there is a hit, coefficients are extracted for the particular species, similarly local temperature range for species is extracted.
7) A for loop, if-else combination is used to calculate the properties in two ranges, Low temperature and high temperature
for i=1:length(Temp_array)
if(Temp_array(i)<Loc_temp_transi)
cp(i) = R*(C(1)+(C(2)*Temp_array(i))+(C(3)*((Temp_array(i))^2)) + (C(4)*((Temp_array(i))^3)) + (C(5)*((Temp_array(i))^4)));
h(i) = (R*Temp_array(i))*(C(1)+(C(2)*Temp_array(i)/2)+(C(3)*((Temp_array(i))^2)/3) + (C(4)*((Temp_array(i))^3)/4) + (C(5)*((Temp_array(i))^4)/5)+(C(6)/Temp_array(i)));
s(i) = R*((C(1)*log(Temp_array(i)))+(C(2)*Temp_array(i))+(C(3)*((Temp_array(i))^2)/2) + (C(4)*((Temp_array(i))^3)/3) + (C(5)*((Temp_array(i))^4)/4)+C(7));
ite = ite +1;
else
cp(i) = R*(C(8)+(C(9)*Temp_array(i))+(C(10)*((Temp_array(i))^2)) + (C(11)*((Temp_array(i))^3)) + (C(12)*((Temp_array(i))^4)));
h(i) = (R*Temp_array(i))*(C(8)+(C(9)*Temp_array(i)/2)+(C(10)*((Temp_array(i))^2)/3) + (C(11)*((Temp_array(i))^3)/4) + (C(12)*((Temp_array(i))^4)/5)+(C(13)/Temp_array(i)));
s(i) = R*((C(8)*log(Temp_array(i)))+(C(9)*Temp_array(i))+(C(10)*((Temp_array(i))^2)/2) + (C(11)*((Temp_array(i))^3)/3) + (C(12)*((Temp_array(i))^4)/4)+C(14));
endif
endfor8) The results are then plotted, images are saved in the directory.
mkdir(['J:\Octave\NASA\',User_species]);
directory = (['J:\Octave\NASA\',User_species]);
cd(directory)
%Plotting
figure(1)
plot(Temp_array(1:ite),cp(1:ite),'color','b',Temp_array(ite:y),cp(ite:y),'color','r')
xlabel('Temperature(K)')
ylabel('Specific heat(kJ/kg-K)')
title('Specific heat')
grid on
legend('Low Temperature','High Temperature')
figure(2)
plot(Temp_array(1:ite),h(1:ite),'color','g',Temp_array(ite:y),h(ite:y),'color','m')
xlabel('Temperature(K)')
ylabel('Enthalpy(kJ)')
title('Enthalpy')
grid on
legend('Low Temperature','High Temperature')
figure(3)
plot(Temp_array(1:ite),s(1:ite),'color','k',Temp_array(ite:y),s(ite:y),'color','c')
xlabel('Temperature(K)')
ylabel('Entropy(kJ/K)')
title('Entropy')
grid on
legend('Low Temperature','High Temperature')
fprintf('Properties calculated for %s\n',B{1})
saveas(1,'Specific_heat.jpg')
saveas(2,'Enthalpy.jpg')
saveas(3,'Entropy.jpg')9) Molecular weight is calculated with the combination of FOR-if loops, linear search is employed
%Molecular weight calculation
mol_wt = mol(User_species);10) The while loop repeats until either the species is exhausted or till there is a hit.
Modifications to move to MODE 2
To move from Mode 1 to Mode 2, user input has to be disabled, if - else condition comparing the Strings has to be disabled, statements in the else section has to be disabled barring the counter.
Loctemp = fgetl(nasa);
B = strsplit(Loctemp,' ');
Loc_species = B{1};
%w = strcmp(User_species,Loc_species);
%if(w==1)
%loop = 0;
%for i=6:(length(B)-1)
%Loc_temp(i-5) = str2num(B{i});
%end
.
.
.
.
..
%else
count = count+1;
loop = 1;
%fgetl(nasa);
%fgetl(nasa);
%fgetl(nasa);
%endifComplete code is as shown below:
%This is coded by Sweekar, [email protected], please contact for queries
%This is an intellectual property of Sweekar, conatct before copying
%This code is used to read from NASA coefficients file for different species
%The code allows user to enter the desired species for which the calculation required
%User is expected to enter the species name by its chemical convention
%Code gets the coefficients from the NASA file and determines enthalpy, entropy and specific heat
%Plots are created for easy visualisation of calculation
%This code will take user inputs, generate plots, calculate molecular weight of species, create the directory and store the plots
%A minor change to this code will enable in generating plots, directories, molecular weight for all the species in Thermo.dat file
%cleanup
close all
clear all
clc
R = 8.314; %Universal gas constant
num_species = 53; %Number of species in the NASA dat file
%Input from user
fprintf("..............NASA Coefficients and Properties calc...............\n")
fprintf("\n The NASA file referenced has %d species\n",num_species)
fprintf("DISCLAIMER: If no plot is obtained, then species is not listed in dat file\n")
fprintf("Please change the details and specify the number of species accurately\n")
fprintf("Plesae enter the desired species name by its chemical formula\n") %Comment this line if program is used to generate properties of all species in dat file
fprintf("Enter within quotes, ex, if oxygen is required enter 'o2'\n")%Comment this line if program is used to generate properties of all species in dat file
User_species = input("Enter species of interest within single quotes: ");%Comment this line if program is used to generate properties of all species in dat file
User_species = upper(User_species); %Making case insensitive
nasa = fopen('THERMO.dat','r'); %Reading Nasa file
fgetl(nasa);
%Global Temperature extraction
A = fgetl(nasa);
Temp = strsplit(A,' ');
for i=2:length(Temp)
Glob_temp(i-1) = str2num(Temp{i});
end
Glob_temp
%Headings extraction
for j = 1:3
fgetl(nasa);
endfor
%variables for Loop and logic execution
loop = 1; %Comment this line if program is used to generate properties of all species in dat file
count = 1;
%To calculated the parameters for the species user desires
while loop == 1 && count<=num_species %change this to count<=53 for second type of program to run
Loctemp = fgetl(nasa);
B = strsplit(Loctemp,' ');
Loc_species = B{1};
w = strcmp(User_species,Loc_species); %Comment this line if program is used to generate properties of all species in dat file
if(w==1) %Comment this line if program is used to generate properties of all species in dat file
loop = 0;
%for i=6:(length(B)-1)
%Loc_temp(i-5) = str2num(B{i});
%end
Loc_temp_low = str2num(Loctemp(49:56))
Loc_temp_high = str2num(Loctemp(58:65))
Loc_temp_transi = str2num(Loctemp(68:75))
Temp_array = [Loc_temp_low:100:Loc_temp_high];
%Coefficient of species extraction
te1 = fgetl(nasa);
te2 = fgetl(nasa);
te3 = fgetl(nasa);
s1 = findstr(te1,'E');
s2 = findstr(te2,'E');
s3 = findstr(te3,'E');
coeffs(1) = str2num(te1(1:s1(1)+3));
for n=2:length(s1)
coeffs(n)=str2num(te1(s1(n-1)+4:s1(n)+3));
end
coeffs1(1) = str2num(te2(1:s2(1)+3));
for n=2:length(s2)
coeffs1(n)=str2num(te2(s2(n-1)+4:s2(n)+3));
end
coeffs2(1) = str2num(te3(1:s3(1)+3));
for n=2:length(s3)
coeffs2(n)=str2num(te3(s3(n-1)+4:s3(n)+3));
end
C = [coeffs coeffs1 coeffs2];
ite = 1;
y = length(Temp_array);
for i=1:length(Temp_array)
if(Temp_array(i)<Loc_temp_transi)
cp(i) = R*(C(1)+(C(2)*Temp_array(i))+(C(3)*((Temp_array(i))^2)) + (C(4)*((Temp_array(i))^3)) + (C(5)*((Temp_array(i))^4)));
h(i) = (R*Temp_array(i))*(C(1)+(C(2)*Temp_array(i)/2)+(C(3)*((Temp_array(i))^2)/3) + (C(4)*((Temp_array(i))^3)/4) + (C(5)*((Temp_array(i))^4)/5)+(C(6)/Temp_array(i)));
s(i) = R*((C(1)*log(Temp_array(i)))+(C(2)*Temp_array(i))+(C(3)*((Temp_array(i))^2)/2) + (C(4)*((Temp_array(i))^3)/3) + (C(5)*((Temp_array(i))^4)/4)+C(7));
ite = ite +1;
else
cp(i) = R*(C(8)+(C(9)*Temp_array(i))+(C(10)*((Temp_array(i))^2)) + (C(11)*((Temp_array(i))^3)) + (C(12)*((Temp_array(i))^4)));
h(i) = (R*Temp_array(i))*(C(8)+(C(9)*Temp_array(i)/2)+(C(10)*((Temp_array(i))^2)/3) + (C(11)*((Temp_array(i))^3)/4) + (C(12)*((Temp_array(i))^4)/5)+(C(13)/Temp_array(i)));
s(i) = R*((C(8)*log(Temp_array(i)))+(C(9)*Temp_array(i))+(C(10)*((Temp_array(i))^2)/2) + (C(11)*((Temp_array(i))^3)/3) + (C(12)*((Temp_array(i))^4)/4)+C(14));
endif
endfor
%Molecular weight calculation
mol_wt = mol(User_species); %Need to change User_species to Loc_species for Mode2
mkdir(['J:\Octave\NASA\',User_species]);
directory = (['J:\Octave\NASA\',User_species]);
cd(directory)
%Plotting
figure(1)
plot(Temp_array(1:ite),cp(1:ite),'color','b',Temp_array(ite:y),cp(ite:y),'color','r')
xlabel('Temperature(K)')
ylabel('Specific heat(kJ/kg-K)')
title('Specific heat')
grid on
legend('Low Temperature','High Temperature')
figure(2)
plot(Temp_array(1:ite),h(1:ite),'color','g',Temp_array(ite:y),h(ite:y),'color','m')
xlabel('Temperature(K)')
ylabel('Enthalpy(kJ)')
title('Enthalpy')
grid on
legend('Low Temperature','High Temperature')
figure(3)
plot(Temp_array(1:ite),s(1:ite),'color','k',Temp_array(ite:y),s(ite:y),'color','c')
xlabel('Temperature(K)')
ylabel('Entropy(kJ/K)')
title('Entropy')
grid on
legend('Low Temperature','High Temperature')
fprintf('Properties calculated for %s\n',B{1})
saveas(1,'Specific_heat.jpg')
saveas(2,'Enthalpy.jpg')
saveas(3,'Entropy.jpg')
cd(['J:\Octave'])
else %Comment this line if program is used to generate properties of all species in dat file
count = count+1;
loop = 1; %Comment this line if program is used to generate properties of all species in dat file
fgetl(nasa); %Comment this line if program is used to generate properties of all species in dat file
fgetl(nasa); %Comment this line if program is used to generate properties of all species in dat file
fgetl(nasa); %Comment this line if program is used to generate properties of all species in dat file
endif %Comment this line if program is used to generate properties of all species in dat file
end
printf('Species number in NASA file %d\n',count)
printf('Species molecular weight is %f \n',mol_wt)
fclose(nasa);
function MW = mol(a)
std_ele = {'O' 'H' 'C' 'N' 'A' 'R'};
std_weights = [16 1 12 14 40 0];
std_numbers = {'1' '2' '3' '4' '5' '6' '7' '8' '9' '10' '20' '30' '80'};
mol_wt = 0;
for i = 1:length(a)
for j = 1:length(std_ele)
if (strcmp(a(i),std_ele{j}))
mol_wt = mol_wt + std_weights(j);
break
endif
endfor
for k = 1:length(std_numbers)
if(strcmp(a(i),std_numbers{k}))
for l = 1:length(std_ele)
if(strcmp(a(i-1),std_ele{l}))
mol_wt = mol_wt + ((str2num(std_numbers{k}))-1)*std_weights(l);
endif
endfor
endif
endfor
endfor
MW = mol_wt;
endfunctionResults:
1) Storing of plots in separate folders in a directory for each species (From Mode 2)
2) Printing of molecular weights of species in the command window
3) Plots for O2,
4) Plots for N2
5) Plots for CO2
Discussion and Report
1) Coefficients are species specific
2) Validation of the results from the plots and the various handbooks return same value
3) The trends of enthalpy, entropy and specific heat for some interesting species such as N2 and CO2 explain concepts of self pressurising gas, crygoneic application
4) Parsing helps read, analyse and utilise processed data very easily
Reference:
https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19940013151.pdf