Parsing NASA thermodynamic data
The main motive of this code is to read the thermodynamic data file to extract the polynomials of the various species and to calculate various properties such as specific heat , enthalpy and entropy.
The formulaes developed by the NASA are given below
where a1,a2,a3,a4,a5,a6,a7 are the coefficients and T are the local temperatures which varies for each species
A small snipet of the thermodynamic data file is shown below
basically the above snippet gives a small idea of how the data looks like and the code is written to read this file line by line.
code:
%project 1
clear all
close all
clc
R=8.31;
%reading the file
f1=fopen(\'thermodata.txt\',\'r\')
line1=fgetl(f1)
%reading the first 4 lines
for i=1:4
fgetl(f1)
end
temp=linspace(200,6000,50)
%reading nthe remaining data
for i=1:52
line=fgetl(f1)
A=strsplit(line ,\' \')
species=A(1)
line2=fgets(f1)
%finding the particular value in the string
sf=strfind(line2,\'E\')
cfh1=line2(1:sf(1)+3);
cfh2=line2(sf(1)+4:sf(2)+3);
cfh3=line2(sf(2)+4:sf(3)+3);
cfh4=line2(sf(3)+4:sf(4)+3);
cfh5=line2(sf(4)+4:sf(5)+3);
line3=fgets(f1)
cfh6=line3(1:sf(1)+3);
cfh7=line3(sf(1)+4:sf(2)+3);
cfL1=line3(sf(2)+4:sf(3)+3);
cfL2=line3(sf(3)+4:sf(4)+3);
cfL3=line3(sf(4)+4:sf(5)+3);
line4=fgets(f1)
cfL4=line4(1:sf(1)+3);
cfL5=line4(sf(1)+4:sf(2)+3);
cfL6=line4(sf(2)+4:sf(3)+3);
cfL7=line4(sf(3)+4:sf(4)+3);
%converting strings to number
hc1=str2num(cfh1);
hc2=str2num(cfh2);
hc3=str2num(cfh3);
hc4=str2num(cfh4);
hc5=str2num(cfh5);
hc6=str2num(cfh6);
hc7=str2num(cfh7);
Lc1=str2num(cfL1);
Lc2=str2num(cfL2);
Lc3=str2num(cfL3);
Lc4=str2num(cfL4);
Lc5=str2num(cfL5);
Lc6=str2num(cfL6);
Lc7=str2num(cfL7);
lc=[Lc1 Lc2 Lc3 Lc4 Lc5 Lc6 Lc7];
hc=[hc1 hc2 hc3 hc4 hc5 hc6 hc7];
for j=1:length(temp)
t=1000;
if temp(j)<=t
Cp(j) = R*(lc(1)+lc(2)*temp(j)+lc(3)*(temp(j))^2 +lc(4)*(temp(j))^3 +lc(5)*(temp(j))^4);
H(j) = R*temp(j)*(lc(1)+lc(2)*temp(j)/2 +lc(3)*((temp(j))^2)/3 +lc(4)*((temp(j))^3)/4 + lc(5)*((temp(5))^4)/5 + lc(6)/temp(j));
S(j) = R*(lc(1)*log(temp(j))+lc(2)*temp(j)+lc(3)*((temp(j))^2)/2 +lc(4)*((temp(j))^3)/3 + lc(5)*((temp(j))^4)/4 + lc(7));
else
Cp(j) = R*(hc(1)+hc(2)*temp(j)+hc(3)*(temp(j))^2 +hc(4)*(temp(j))^3 +hc(5)*(temp(j))^4);
H(j) = R*temp(j)*(hc(1)+hc(2)*temp(j)/2 +hc(3)*((temp(j))^2)/3 +hc(4)*((temp(j))^3)/4 + hc(5)*((temp(5))^4)/5 + hc(6)/temp(j));
S(j) = R*(hc(1)*log(temp(j))+hc(2)*temp(j)+hc(3)*((temp(j))^2)/2 +hc(4)*((temp(j))^3)/3 + hc(5)*((temp(j))^4)/4 + hc(7));
end
end
spec=char(species)
fd = sprintf(\'%s\',spec);
cur_dir=pwd
mkdir(spec);
cd(spec);
%plotting
%Plot 1, Specific heat vs temp
figure(1);
plot(temp,Cp);
xlabel(\'Temperature in K\');
ylabel(\'Specific heat in kJ\')
title(\'Specific heat vs Temp\');
%Plot 2, Enthalpy vs temp
figure(2);
plot(temp,H);
xlabel(\'Temperature in K\');
ylabel(\'Enthalpy in kJ/(mol)\')
title(\'Enthalpy vs Temp\');
%Plot 3, Entropy vs temp
figure(3);
plot(temp,S);
xlabel(\'Temperature in K\');
ylabel(\'Entropy in kJ/(mol K)\')
title(\'Entropy vs Temp\');
% saving the plots
saveas(1, \'cp.png\');
saveas(2, \'h.png\');
saveas(3, \'s.png\');
cd(cur_dir);
end
cfh=coefficient for higher temperature range(first 7)
cfL=coefficient of lower temperature range(last 7)
individual files are made in the same folder by using the mkdir command which means (make directory)
cd(specs) means making the individual files in the current directory.
rest of the code is self explanatory
screenshot:
ouput of AR:
1.Temp vs Cp
2.Temp vs Enthalpy(H)
3.Temp vs Entropy(S)
link to access all the folder is given below
https://drive.google.com/drive/folders/1gSa7ITxIzGKjoYjY4_Bodq5f101HwfKQ?usp=sharing