Project 1 - Parsing NASA thermodynamic data
Aim
To write a function that extracts the 14 co-efficients and to calculate the enthalpy, entropy and specific heats for all the species in the data file.
Formulae Used
`Cp=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*lnT+a2*T+a3*T^2/2+a4*T^3/3+a5*T^4/4+a7`
where Cp is specific heat
H is enthalpy
S is entropy
R is gas constant
T is temperature
The first 7 coefficients are high-temperature coefficients and the second 7 coefficients are low-temperature coefficients.
Procedure:
- First we have to find the global temperatures from the data file)
- Then we have to extract the coefficients from the file for all the species
- Using the formulae, calculate the Entropy, enthalpy and specific heat of the species
- Plot the results in the respective folder name of the species
- Determining the molecular weight of the species
Coding and explanation:
clear all
close all
clc
fileID = fopen('exp.txt','w') To create a new text document to store the molecular masses
f1=fopen('THERMO.dat','r') Opens the thermo data file
fgetl(f1) executes the next line
tline1=fgetl(f1) This line contains global temperatures. we splits the line in to cell arrays
tline2=strsplit(tline1,'') by splitting the values seperated by space.
temp=str2num(tline2{1}) Cell arrays are converted into array consisting of numbers
global_low_temp=temp(1,1)
global_mid_temp=temp(1,2) Global low mid & high temperatures are extracted from the arrays
global_high_temp=temp(1,3)
tline3=fgetl(f1);
tline4=fgetl(f1);
tline5=fgetl(f1);
for i=1:53 Using for loop, we run the loop for 53 number of species
tline6=fgetl(f1) O L 1/90O 1 G 200.000 3500.000 1000.000 1
x=findstr(tline6,'G') Here the O is the species name, the L 1 /900 is not necessary the
local_temp=tline6(x:end) remaining is the local temperatures. We extract them by storing in cell
local_temp=strsplit(local_temp) arrays.
local_low_temp=local_temp(1,2); Similarly we convert the cell arrays into numbers
local_low_temp=str2num(local_low_temp{1})
local_mid_temp=local_temp(1,4);
local_mid_temp=str2num(local_mid_temp{1})
local_high_temp=local_temp(1,3);
local_high_temp=str2num(local_high_temp{1})
2.56942078E+00-8.59741137E-05 4.19484589E-08-1.00177799E-11 1.22833691E-15 2 tline7=fgetl(f1)
a=findstr(tline7,'E'); Here we find the places of E , the digits which follows E is the power of 10.
a1=tline7(1:a(1)+3); after E there are 3 places which represents the power of 10. we store them in cell
a1=str2num(a1) and convert them in to numbers
a2=tline7(a(1)+4:a(2)+3);
a2=str2num(a2) Using the same logic we extract all the 14 coefficients
a3=tline7(a(2)+4:a(3)+3);
a3=str2num(a3)
a4=tline7(a(3)+4:a(4)+3);
a4=str2num(a4)
a5=tline7(a(4)+4:a(5)+3);
a5=str2num(a5)
tline8=fgetl(f1);
b=findstr(tline8,'E');
a6=tline8(1:b(1)+3);
a6=str2num(a6)
a7=tline8(b(1)+4:b(2)+3);
a7=str2num(a7)
a8=tline8(b(2)+4:b(3)+3);
a8=str2num(a8)
a9=tline8(b(3)+4:b(4)+3);
a9=str2num(a9)
a10=tline8(b(4)+4:b(5)+3);
a10=str2num(a10)
tline9=fgetl(f1)
c=findstr(tline9,'E');
a11=tline9(1:c(1)+3);
a11=str2num(a11)
a12=tline9(c(1)+4:c(2)+3);
a12=str2num(a12)
a13=tline9(c(2)+4:c(3)+3);
a13=str2num(a13)
a14=tline9(c(3)+4:c(4)+3);
a14=str2num(a14)
% end
R=8.314; R is universal gas constant
T=local_low_temp:local_high_temp;
T1=local_low_temp
T2=local_mid_temp
T3=local_high_temp
`Cp=a1+a2*T+a3*T^2+a4*T^3+a5*T^4`
if T<T2 for low temperature , we use coefficients a1 to a7
Cp=(a1+a2*T+a3*T.^2+a4*T.^3+a5*T.^4)*R
else T>T2 for high temperature , we use coefficients a8 to a14
Cp=(a7+a8*T+a9*T.^2+a10*T.^3+a11*T.^4)*R
end
if T<T2 `H/(R*T)=a1+a2*T/2+a3*T^2/3+a4*T^3/4+a5*T^4/5+(a6)/T`
H=R*T.*(a1+(a2/2)*T+(a3/3)*T.^2+(a4/4)*T.^3+(a5/5)*T.^4+a6*(1.\T));
else T>T2
H=R*T.*(a8+(a9/2)*T+(a10/3)*T.^2+(a11/4)*T.^3+(a12/5)*T.^4+a13*(1.\T));
end
`S/R=a1*lnT+a2*T+a3*T^2/2+a4*T^3/3+a5*T^4/4+a7`
if T<T2
S=R*(a1*log(T)+a2*T+(a3/2)*T.^2+(a4/3)*T.^3+(a5/4)*T.^4+a7);
else T>T2
S=R*(a8*log(T)+a9*T+(a10/2)*T.^2+(a11/3)*T.^3+(a12/4)*T.^4+a14);
end
figure(1) plotting the Temperature vs. Specific heat
plot(T,Cp)
ylabel('Cp')
xlabel('Temperature')
figure(2) plotting the Temperature vs. Enthalpy
plot(T,H)
ylabel('Enthalpy')
xlabel('Temperature')
figure(3) plotting the Temperature vs. Entropy
plot(T,S)
ylabel('Entropy')
xlabel('Temperature')
spec_name=strsplit(tline6) Extracting the species name
spec_name=spec_name{1}
mkdir(['F:\siva raam\Mat Prac\File_parsing\plots\',spec_name]) mkdir is used to create a new folder directory
cd(['F:\siva raam\Mat Prac\File_parsing\plots\',spec_name]) cd is used to change the directory
saveas(1,'specific heat','jpg')
saveas(2,'enthalpy','jpg')
saveas(3,'entropy','jpg')
atom = ['H','C','O','N','A']; we know the species are having elements H,C,O,N,A
atm_weight = [1.00784 12.0107 15.999 14.0067 39.948]; and the individual atomic weight are known
M_Weight = 0;
for i = 1:length(spec_name) we use two for loops, first loop for number of elements
for j = 1:length(atom) in the species & the second loop runs for total number of elements
if strcmp(spec_name(i),atom(j))
M_Weight = M_Weight +atm_weight(j); store values if character match.
position = j ;%store position
end
end
n = str2double(spec_name(i)); for multiplication of the repeating elements
if n>1
M_Weight = M_Weight + atm_weight(position)*(n-1);
end
%display molecular weight
fprintf(fileID,'Molecular Weight of %s : %f',spec_name,M_Weight) printing the data nad saving in a text file
fprintf(fileID,'\n')
end
end
clear all
close all
clc
fileID = fopen('exp.txt','w')
f1=fopen('THERMO.dat','r')
fgetl(f1)
tline1=fgetl(f1)
tline2=strsplit(tline1,'')
temp=str2num(tline2{1})
global_low_temp=temp(1,1)
global_mid_temp=temp(1,2)
global_high_temp=temp(1,3)
tline3=fgetl(f1);
tline4=fgetl(f1);
tline5=fgetl(f1);
for i=1:53
tline6=fgetl(f1)
x=findstr(tline6,'G')
local_temp=tline6(x:end)
local_temp=strsplit(local_temp)
local_low_temp=local_temp(1,2);
local_low_temp=str2num(local_low_temp{1})
local_mid_temp=local_temp(1,4);
local_mid_temp=str2num(local_mid_temp{1})
local_high_temp=local_temp(1,3);
local_high_temp=str2num(local_high_temp{1})
tline7=fgetl(f1)
a=findstr(tline7,'E');
a1=tline7(1:a(1)+3);
a1=str2num(a1)
a2=tline7(a(1)+4:a(2)+3);
a2=str2num(a2)
a3=tline7(a(2)+4:a(3)+3);
a3=str2num(a3)
a4=tline7(a(3)+4:a(4)+3);
a4=str2num(a4)
a5=tline7(a(4)+4:a(5)+3);
a5=str2num(a5)
tline8=fgetl(f1);
b=findstr(tline8,'E');
a6=tline8(1:b(1)+3);
a6=str2num(a6)
a7=tline8(b(1)+4:b(2)+3);
a7=str2num(a7)
a8=tline8(b(2)+4:b(3)+3);
a8=str2num(a8)
a9=tline8(b(3)+4:b(4)+3);
a9=str2num(a9)
a10=tline8(b(4)+4:b(5)+3);
a10=str2num(a10)
tline9=fgetl(f1)
c=findstr(tline9,'E');
a11=tline9(1:c(1)+3);
a11=str2num(a11)
a12=tline9(c(1)+4:c(2)+3);
a12=str2num(a12)
a13=tline9(c(2)+4:c(3)+3);
a13=str2num(a13)
a14=tline9(c(3)+4:c(4)+3);
a14=str2num(a14)
% end
R=8.314;
T=local_low_temp:local_high_temp;
T1=local_low_temp
T2=local_mid_temp
T3=local_high_temp
if T<T2
Cp=(a1+a2*T+a3*T.^2+a4*T.^3+a5*T.^4)*R
else T>T2
Cp=(a7+a8*T+a9*T.^2+a10*T.^3+a11*T.^4)*R
end
if T<T2
H=R*T.*(a1+(a2/2)*T+(a3/3)*T.^2+(a4/4)*T.^3+(a5/5)*T.^4+a6*(1.\T));
else T>T2
H=R*T.*(a8+(a9/2)*T+(a10/3)*T.^2+(a11/4)*T.^3+(a12/5)*T.^4+a13*(1.\T));
end
if T<T2
S=R*(a1*log(T)+a2*T+(a3/2)*T.^2+(a4/3)*T.^3+(a5/4)*T.^4+a7);
else T>T2
S=R*(a8*log(T)+a9*T+(a10/2)*T.^2+(a11/3)*T.^3+(a12/4)*T.^4+a14);
end
figure(1)
plot(T,Cp)
ylabel('Cp')
xlabel('Temperature')
figure(2)
plot(T,H)
ylabel('Enthalpy')
xlabel('Temperature')
figure(3)
plot(T,S)
ylabel('Entropy')
xlabel('Temperature')
spec_name=strsplit(tline6)
spec_name=spec_name{1}
mkdir(['F:\siva raam\Mat Prac\File_parsing\plots\',spec_name])
cd(['F:\siva raam\Mat Prac\File_parsing\plots\',spec_name])
saveas(1,'specific heat','jpg')
saveas(2,'enthalpy','jpg')
saveas(3,'entropy','jpg')
atom = ['H','C','O','N','A'];
atm_weight = [1.00784 12.0107 15.999 14.0067 39.948];
M_Weight = 0;
% current_species =char(spec_name);
for i = 1:length(spec_name)
for j = 1:length(atom)
if strcmp(spec_name(i),atom(j))
M_Weight = M_Weight +atm_weight(j);%store values if character match.
position = j ;%store position
end
end
n = str2double(spec_name(i));
if n>1
M_Weight = M_Weight + atm_weight(position)*(n-1);
end
%display molecular weight
fprintf(fileID,'Molecular Weight of %s : %f',spec_name,M_Weight)
fprintf(fileID,'\n')
end
end
Result:
The plots temperature Cp, H & S are plotted for all the species and the corresponding molecular masses are determined