Parsing NASA thermodynamic data using MATLAB
Aim
The aim of this program is to parse the given data file which contains thermodynamic coefficients of various species and calculate their specific heat, enthalpy and entropy values over a temperature range.
Formulae used
Equations used in this program are:
- `Cp=R*(a1+a2*T+a3*T^2+a4*T^3+a5*T^4)`
- `H=R∗T∗(a1+a2∗T/2+a3∗T^2/3+a4∗T^3/4+a5∗T^4/5+(a6)/T)`
- `S=R*( a1*ln(T) +a2*T + a3*T^2/2+ a4*T^3/3 +a5*T^4/4 +a7)`
Where the coefficients are obtained from the given file along with the temperature range for each species.
Program
The file is opened in read mode using 'fopen' command and the fileID is stored in 'file'. The first five lines of the data file are skipped as they are not required. 'fgetl' command is used. Each species has a header line and 3 coefficients lines. If the current line is a header line, the 'Nasacompandtemp_func' is called and if it is a coefficients line the 'NASAcoeffs_func' is called. This is implemented using a while loop and if statement. The species names are stored in a string array 'comp' and the local temperature ranges in a cell array 'temp'.
while(a~="END")
if(mod(index,4)==1) %Parsing header line for each species
[comp(index1),temp{index1}]=NASAcompandtemp_func(a);
index1=index1+1;
else %Parsing the coefficients block for each species
coeffs{index2}=NASAcoeffs_func(a);
index2=index2+1;
end
a=fgetl(file);
index=index+1;
endThe 'NASAcompandtemp_func' splits the header line using 'strsplit' command with ' ' (space) as the delimiter. The species name and the local temperature range are obtained according to their index.
% parses the header line for each species
function [comp,temp]=NASAcompandtemp_func(a)
a=strsplit(a,' ');
comp=string(a(1)); %stores the species name
c=1;
for i=(length(a)-3):(length(a)-1)
temp(c)=str2num(a{i}); %stores the local temperature range of the species
c=c+1;
endThe 'NASAcoeffs_func' uses the 'strfind' command to locate the position of 'E' and splits the current line accordingly to obtain the coefficients and stores them in 'coeffs' which is a cell array.
% parses the coefficients block
function coeffs=NASAcoeffs_func(b)
f=strfind(b,"E");
initpos=1;
for i=1:length(f)
coeffs(i)=str2num(b(initpos:f(i)+3)); %stores the coefficients for each species
initpos=f(i)+4;
endThe molecular weight of each species is calculated through the 'NASAmolwtcal_func' and stored in array 'molwt'. It is simultaneously displayed in the command window. The function uses the logic that if the species is an atom or CH2(S) the atomic weight is directly stored and if it’s a molecule the weight of each atom is added to obtain the molecular weight.
%calculates the molecular weight of each species
function molwt=NASAmolwtcal_func(var)
atoms=["H","C","O","N","AR","CH2(S)"];
atwt=[1,12,16,14,40,14];
molwt=0;
if(ismember(var,atoms))
index=find(atoms==var);
molwt=atwt(index);
else
for i=1:length(var)
if(ismember(var(i),atoms))% if its an atom directly add
index=find(atoms==var(i));
molwt=molwt+atwt(index);
else %if its the number(n) following an atom, add the weight of the previous atom n-1 times
molwt=molwt+atwt(index)*(str2num(var(i))-1);
end
end
end The user is then asked to input the species for which s/he wants the thermodynamic plots for. The user is requested to try another species If the the program doesn't have the data for inputted species. The specific heat, enthalpy and entropy are calculated for the local temperature in the 'NASAcalculate_func'. Once the data is obtained they are plotted and the figure is saved using the 'NASAplotting_func'. The user is then if s/he wants the plots for another species. Based on the users input the program continues or end.
b=1;
b1=1;
% To get input from user and providing the plots for the species requested
while(b1==1)
while(b==1)
species=input("enter name of species\n");
species=upper(species);
if(ismember(species,comp)) %if the program has the data fro inputted species
index=find(comp==species);
[cp,h,s,t]=NASAcalculate_func(index,temp{index},coeffs);
b=NASAplotting_func(index,temp,comp,cp,h,s,t);
else
disp('Try another species');
end
end
resp=input("Enter 'yes' if you want the thermodynamic data for another species\n");
resp=lower(resp);
if(strcmp(resp,'yes'))
b=1;
else
disp('Thanks for using this program')
b1=0;
end
end The 'NASAcalculate_func' calculates the specific heat, enthalpy and entropy values using the equations mentioned above.
%calculates the values of specific heat, enthalpy and entropy for given temperature range
function [cp,h,s,t]=NASAcalculate_func(index,tempvals,coeffs)
lt=tempvals(1);
mt=tempvals(3);
ht=tempvals(2);
c1_5=coeffs{index*3-2};
c6_10=coeffs{index*3-1};
c11_14=coeffs{index*3};
t=linspace(lt,ht,500);
for i=1:500
if(t(i)>mt)
cp(i)=8.314*(c1_5(1) + c1_5(2)*t(i) + c1_5(3)*t(i)^2 + c1_5(4)*t(i)^3 + c1_5(5)*t(i)^4);
h(i)= 8.314*t(i)*(c1_5(1) + c1_5(2)*t(i)/2 + c1_5(3)*t(i)^2/3 + c1_5(4)*t(i)^3/4 + c1_5(5)*t(i)^4/5 + c6_10(1)/t(i));
s(i)= 8.314*(c1_5(1)*log(t(i)) + c1_5(2)*t(i) + c1_5(3)*t(i)^2/2 + c1_5(4)*t(i)^3/3 + c1_5(5)*t(i)^4/4 + c6_10(2));
else
cp(i)=8.314*(c6_10(3) + c6_10(4)*t(i) + c6_10(5)*t(i)^2 + c11_14(1)*t(i)^3 + c11_14(2)*t(i)^4);
h(i)= 8.314*t(i)*(c6_10(3) + c6_10(4)*t(i)/2 + c6_10(5)*t(i)^2/3 + c11_14(1)*t(i)^3/4 + c11_14(2)*t(i)^4/5 + c11_14(3)/t(i));
s(i)= 8.314*(c6_10(3)*log(t(i)) + c6_10(4)*t(i) + c6_10(5)*t(i)^2/2 + c11_14(1)*t(i)^3/3 + c11_14(2)*t(i)^4/4 + c11_14(4));
end
endThe NASAplotting_func uses the 'mkdir' command to create a new folder in the current directory and the 'saveas' command to save the plots for the requested species.
% plots the specific heat, enthalpy an dentropy against temperature for each species
%and stores the figure as png file
function b=plotting_func(index,temp,comp,cp,h,s,t)
mkdir NASA_dataparsing_images;
pathname="C:\Users\Venkat JR\Documents\MATLAB\NASA_dataparsing_images";
mkdir (pathname,comp{index});
newpathname=pathname+"\"+comp{index};
figure(1)
plot(t,cp,'linewidth',3);
xlim([0,temp{index}(2)+200]);
xlabel('Temperature(K)');
ylabel('Specific Heat capacity(J/kg-K)')
title(['Species: ',comp{index}]);
saveas(1, fullfile(newpathname,'CPvsT'), 'png');
figure(2)
plot(t,h,'linewidth',3);
xlim([0,temp{index}(2)+200]);
xlabel('Temperature(K)');
ylabel('Enthalpy(J/kg)');
title(['Species: ',comp{index}]);
saveas(2, fullfile(newpathname,'HvsT'), 'png');
figure(3)
plot(t,s,'linewidth',3);
xlim([0,temp{index}(2)+200]);
xlabel('Temperature(K)');
ylabel('Entropy(J/K)');
title(['Species: ',comp{index}]);
saveas(3, fullfile(newpathname,'SvsT'), 'png');
b=0;
endThe following bit of code is used to create plots for every species and store them in separate folder. It is commented out in the program so it doesn’t increase the running time.
% The following block was used to save the plots for every species and store them in separate folders
%{
for i=1:length(comp)
[cp,h,s,t]=calculate_func(i,temp{i},coeffs);
b=NASAplotting_func(i,temp,comp,cp,h,s,t);
end
%}The file is closed using 'fclose' command.
While running the program a warning might occur if the folder being created in the 'NASAplotting_func' to store the plot figures, already exists in the current directory.
Main Code
clear all
close all
clc
file=fopen('THERMO.dat','r');
for i=1:5 %Skipping the first five lines
fgetl(file);
end
a=fgetl(file);
index=1;
index1=1;
index2=1;
while(a~="END")
if(mod(index,4)==1) %Parsing header line for each species
[comp(index1),temp{index1}]=NASAcompandtemp_func(a);
index1=index1+1;
else %Parsing the coefficients block for each species
coeffs{index2}=NASAcoeffs_func(a);
index2=index2+1;
end
a=fgetl(file);
index=index+1;
end
disp('Molecular weight of every species is displayed below');
for i=1:length(comp) % Calculating the molecular wieght for each species and displaying them
molwt(i)=NASAmolwtcal_func(comp{i});
disp(comp{i});
disp(molwt(i));
end
% The following block was used to save the plots for every species and store them in separate folders
%{
for i=1:length(comp)
[cp,h,s,t]=calculate_func(i,temp{i},coeffs);
b=NASAplotting_func(i,temp,comp,cp,h,s,t);
end
%}
b=1;
b1=1;
% To get input from user and providing the plots for the species requested
while(b1==1)
while(b==1)
species=input("enter name of species\n");
species=upper(species);
if(ismember(species,comp)) %if the program has the data fro inputted species
index=find(comp==species);
[cp,h,s,t]=NASAcalculate_func(index,temp{index},coeffs);
b=NASAplotting_func(index,temp,comp,cp,h,s,t);
else
disp('Try another species');
end
end
resp=input("Enter 'yes' if you want the thermodynamic data for another species\n");
resp=lower(resp);
if(strcmp(resp,'yes'))
b=1;
else
disp('Thanks for using this program')
b1=0;
end
end
fclose(file);Results
Plots for O2
Plots for N2
Plots for CO2
References