Parsing NASA thermodynamic data
AIM:
Parsing the NASA THERMODYNAMIC DATA.
OBJECTIVES OF THE PROJECT:
A function that extracts the 14 coefficients and calculates the enthalpy, entropy and specific heats for all the species in the data file.
Calculate the molecular weight of each species and display it in the command window.
Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature: high temperature) specific for each species.
FORMULA USED:
`C_P/R=a_1+a_2 T+a_3 T^2+a_4 T^3+a_5 T^4`
`H/(RT)=a_1+(a_2 T)/2+(a_3 T^2)/3+(a_4 T^3)/4+(a_5 T^4)/5+a_6/T`
`s/R=a_1 ln(T)+a_2 T+(a_3 T^2)/2+(a_4 T^3)/3+(a_5 T^4)/4+a_7`
INTRODUCTION:
- Effect of temperature on specific heat
Specific heat is a measure of the ability of the substance to absorb heat.
The heat goes first into increasing the kinetic energies of the molecules.
Molecules can also store energy in vibrations and rotations. These energies are quantized
For example, the O-H bonds in water can stretch and bend, and the whole molecule can rotate about its x-, y-, and z-axes.
At low temperatures, collisions do not provide enough energy to get out of the ground states for rotation or vibration.
The average energy increases only from translation.
As the substance heats up, the average kinetic energy of the molecules increases.
The collisions impart enough energy to allow rotation to occur.
Rotation then contributes to the internal energy and raises the specific heat.
The vibrational energy states are farther apart than the rotation energy states.
Vibration does not jump out of the ground state until high temperatures (≈ 1000 K).
It then contributes to the specific heat.
- Effect of temperature on Enthalpy
When a process occurs at constant pressure, the heat evolved (either released or absorbed) is equal to the change in enthalpy. Enthalpy (H">H) is the sum of the internal energy (U">U) and the product of pressure and volume (PV">PV) given by the equation:
`H=U+PV`.....(1)
When a process occurs at constant pressure, the heat evolved (either released or absorbed) is equal to the change in enthalpy. Enthalpy is a state function which depends entirely on the state functions T, P and U. Enthalpy is usually expressed as the change in enthalpy (ΔH) for a process between initial and final states:
ΔH=ΔU+ΔPV.....(2)
If temperature and pressure remain constant through the process and the work is limited to pressure-volume work, then the enthalpy change is given by the equation:
`ΔH=ΔU+PΔV`...(3)
Also at constant pressure the heat flow (q">q) for the process is equal to the change in enthalpy defined by the equation:
`ΔH=q`.....(4)
When the temperature increases, the amount of molecular interactions also increases. When the number of interactions increase, then the internal energy of the system rises. According to the first equation given, if the internal energy (U">U) increases then the ΔH">ΔH increases as temperature rises. We can use the equation for heat capacity and Equation 2 to derive this relationship.
`C=q/(ΔT)`.....(5)
Under constant pressure, substitute Equation 6">4 into Equation 7">5:
7">`c_p=((ΔH)/(ΔT))_p`....(6)
where the subscript P">P indicates the derivative is done under constant pressure.
- Effect of temperature on Entropy:
Entropy is a measure of the random activity in a system. The entropy of a system depends on your observations at one moment.
When we say random, we mean energy that can't be used for any work.
The formula (delta)S = (delta)Q /(delta)T. "S" is the entropy value, "Q" is the measure of heat, and "T" is the temperature of the system measured in Kelvin degrees.
Several factors affect the amount of entropy in a system. If you increase temperature, you increase entropy.
(1) More energy put into a system excites the molecules and the amount of random activity.
(2) As a gas expands in a system, entropy increases. This one is also easy to visualize. If an atom has more space to bounce around, it will bounce more. Gases and plasmas have large amounts of entropy when compared to liquids and solids.
(3) When a solid becomes a liquid, its entropy increases.
(4) When a liquid becomes a gas, its entropy increases. We just talked about this idea. If you give atoms more room to move around, they will move. You can also think about it in terms of energy put into a system. If you add energy to a solid, it can become a liquid. Liquids have more energy and entropy than solids.
(5) Any chemical reaction that increases the number of gas molecules also increases entropy. A chemical reaction that increases the number of gas molecules would be a reaction that pours energy into a system. More energy gives you greater entropy and randomness of the atoms.
BODY OF THE CONTENT:
The program is basically divided into 5 major parts.
- Reading the Header information:
close all
clear
clc
%% Opening The 'THERMO.dat' File
f1=fopen('THERMO.dat','r');
% Universal gas constant
R=8.314;
%% To Read The Header Information
for i=1:5
if i==2
tempratures=fgetl(f1);
temp=strsplit(tempratures,'');
Tmp=str2num(temp{1});
lower_global_tmp=Tmp(1);
mid_gloabal_tmp=Tmp(2);
higher_gloabal_tmp=Tmp(3);
else
comments{i}=fgetl(f1());
end
endAs the file ‘THERMO.dat’ contains 5 lines of comments in which the 2nd line is the global temperatures.
Firstly opening the file through ‘fopen’ command and than reading the 1st header info and saving the global temperatures(in form of numbers) as well as the comments (in a cell array).
The temperatures are splitted using the command ‘strsplit’. And then converting it into a number using command ‘str2num’ and storing each value into a particular variable.
(Note: while converting the cell array to number array curly braces{} are used as cell is a special type of array which can hold both number as well as string, therefore {} are used to get the access inside the content of the cell array)
- Storing the remaining data into a cell array and dividing them accordance with the need
%% Seprating The Sepcies coloumn and the coefficients coloumn
coefficients=A(:,2:4); % coefficeints columns
num_values=coefficients'; %transposing the cell so that each species coffiecient
%appears in a single row
species=A(:,1); % getting the species column
%converting it into char array
sp=char(species);
%% evaluating the local temperatures as well as creating dividing the rows of each species
for m=1:53
% getting the local tempratures
values{m}=strsplit(sp(m,1:end-1),'G ');
L_tmp=values{m};
local_tmp(m)=L_tmp(1,2);
%for creating species
elements_split{m}=strsplit(sp(m,:),' ');
each_species=elements_split{m};
Species(m)=each_species(1);
%seprating coefficients of each line for every species
w=num_values(:,m);
line1_species(m)=w(1,1);
line2_species(m)=w(2,1);
line3_species(m)=w(3,1);
end
L_Temp=char(local_tmp);
local_tempreatures=str2num(L_Temp);
lower_local_temp=local_tempreatures(:,1);
higher_local_temp=local_tempreatures(:,2);
mid_local_temp=local_tempreatures(:,3);
As after the 5th line each species has its own 4 lines and there are 53 species in total, therefore utilizing the double for loop which will store the info of every species in a set of rows as shown below.
Form the above figure it is seen that each row represents the information of every different species.
Now saving the coefficients and species in a different array’s and utilizing these values to evaluate the name of the species, local temperatures as well as the coefficients values.
(Note: for convenience, the cell array is converted into a char array using the command ‘char’.)
In this 1st splitting the species array in such a way to extract the local temperatures and then splitting it to get the name of the species and lastly separating each line of every species into three different variables.
With the use of the local_tmp variable extracting the temperatures in line 67 to 71.
- Extracting the coefficients from every line generated
%% extracting cofficents from lines of every species
%converting it into character form so that we can use
%the function strfind
line1=char(line1_species);
line2=char(line2_species);
line3=char(line3_species);
for s=1:53
a=strfind(line1(s,:),'E');
% upper tmperature cofficients
Ah1(s,:)=line1(s,1:a(1)+3);
ah1(s)=str2double(Ah1(s,:));
Ah2(s,:)=line1(s,a(1)+4:a(2)+3);
ah2(s)=str2double(Ah2(s,:));
Ah3(s,:)=line1(s,a(2)+4:a(3)+3);
ah3(s)=str2double(Ah3(s,:));
Ah4(s,:)=line1(s,a(3)+4:a(4)+3);
ah4(s)=str2double(Ah4(s,:));
Ah5(s,:)=line1(s,a(4)+4:a(5)+3);
ah5(s)=str2double(Ah5(s,:));
Ah6(s,:)=line2(s,1:a(1)+3);
ah6(s)=str2double(Ah6(s,:));
Ah7(s,:)=line2(s,a(1)+4:a(2)+3);
ah7(s)=str2double(Ah7(s,:));
%lower temperature coffiecients
AL1(s,:)=line2(s,a(2)+4:a(3)+3);
aL1(s)=str2double(AL1(s,:));
AL2(s,:)=line2(s,a(3)+4:a(4)+3);
aL2(s)=str2double(AL2(s,:));
AL3(s,:)=line2(s,a(4)+4:a(5)+3);
aL3(s)=str2double(AL3(s,:));
AL4(s,:)=line3(s,1:a(1)+3);
aL4(s)=str2double(AL4(s,:));
AL5(s,:)=line3(s,a(1)+4:a(2)+3);
aL5(s)=str2double(AL5(s,:));
AL6(s,:)=line3(s,a(2)+4:a(3)+3);
aL6(s)=str2double(AL6(s,:));
AL7(s,:)=line3(s,a(3)+4:a(4)+3);
aL7(s)=str2double(AL7(s,:));
% aligning the contents of each species
c(s,:)=[ah1(s) ah2(s) ah3(s) ah4(s) ah5(s) ah6(s) ah7(s) aL1(s) aL2(s) aL3(s) aL4(s) aL5(s) aL6(s) aL7(s)];
endFirstly, converting the cell array into the character array.
As every line has character ‘E’ in it. Using it as a flag to access the right and left hand side of the string.
(Note: ‘strfind’ command is used to find the character inside a string and it can only find the value one row at a time, value of a is also same for each and every line therefore utilizing the value of a of line 1 only for every other line)
Therefore, once the position of ‘E’ is been searched using it to split the string and converting them into numbers using str2dobule (for more precision) and at each iteration the coefficients of one species is been calculated.
At last aligning the coefficients in a way such that every row represents all 14 coefficients of each species as shown below.
- Calculating the values of Cp, H and S and plotting graph and saving it into a folder named by species.
for p=1:53
T=lower_local_temp(p):higher_local_temp(p);
if T<=mid_local_temp(p)
Cp=R*(aL1(p)+(aL2(p)*T)+(aL3(p)*(T.^2))+(aL4(p)*(T.^3))+(aL5(p)*(T.^4)));
H=R*((aL1(p)*T)+((aL2(p)*(T.^2))/2)+((aL3(p)*(T.^3))/3)+((aL4(p)*(T.^4))/4)+((aL5(p)*(T.^5))/5)+aL6(p));
S=R*((aL1(p)*log(T)) + (aL2(p)*T) + ((aL3(p)*(T.^2))/2) + ((aL4(p)*(T.^3))/3)+((aL5(p)*(T.^4))/4)+(aL7(p)));
else
Cp=R*(ah1(p)+(ah2(p)*T)+(ah3(p)*(T.^2))+(ah4(p)*(T.^3))+(ah5(p)*(T.^4)));
H=R*((ah1(p)*T)+((ah2(p)*(T.^2))/2)+((ah3(p)*(T.^3))/3)+((ah4(p)*(T.^4))/4)+((ah5(p)*(T.^5))/5)+ah6(p));
S=R*((ah1(p)*log(T)) + (ah2(p)*T) + ((ah3(p)*(T.^2))/2) + ((ah4(p)*(T.^3))/3)+((ah5(p)*(T.^4))/4)+(ah7(p)));
end
figure(p)
plot(T,Cp,'r-','linewidth',2);
title(strcat('Specific Heat Vs Temperature of:',Species{:,p}))
xlabel('Temperature(K)')
ylabel('specific Heat(J/mol.K)')
grid on
foldername=sprintf('%s',Species{:,p});
mkdir(foldername(1,:))
picturename = fullfile(foldername(1,:), 'Specific_Heat');
saveas (figure(p), picturename, 'png')
figure(p+53)
plot(T,H,'b-','linewidth',2);
title(strcat('Enthalpy Vs Temperature of:',Species{:,p}))
xlabel('Temperature(K)')
ylabel('Enthalpy(J)')
grid on
picturename = fullfile(foldername(1,:), 'Enthaply');
saveas (figure(p+53), picturename, 'png')
figure(p+106)
plot(T,S,'g-','linewidth',2);
title(strcat('Entropy Vs Temperature of:',Species{:,p}))
xlabel('Temperature(K)')
ylabel('Entropy(J/K)')
grid on
picturename = fullfile(foldername(1,:), 'Entropy');
saveas (figure(p+106), picturename, 'png')
end
As using the formulas for Cp, H and S evaluating its values for the local range of temperatures with the given conditions (1st seven coefficients for higher temp and last 7 coefficients for lower temp).
And then using the basic plot command, the main aim was to save the three plots of each species with the folder named by that species.
With the use of sprint command and the Species array assigning the name of folder and then for making the folder using the command ‘mkdir’ and ‘fullfile’ command is used to direct the picture into the folder created as well as to give name to the picture.
At last using the command ‘saveas’ to save the figure generated into .png form. Same procedure is being applied to all other three graphs.
The result generated is shown below:
- Calculating the value of molecular weight of each species
%% calculation of molecular weight
aw=[16.00,1.00,12.00,14.00,40.00];
elements= ["O","H","C","N","AR"];
for i=1:53
z=char(Species(i));
weight=0;
for j=1:length(z)
for k=1:length(elements)
if strcmpi(z(j),elements(k))
weight=weight+aw(k);
r=k;
else
if str2num(z(j))>=1
val=str2num(z(j));
weight=weight+aw(r)+(aw(r)*(val-2));
break;
end
end
end
end
%specially for AR
if weight==0
weight=aw(k);
end
wt(i)=weight;
disp(strcat('Molecular weight of:',Species{i},'=',num2str(wt(i))));
end(Note: lots of trial and error method have been used in this case to get to the result which was desired.)
To explain the concept used over here, let Us take an example first.
Suppose we need to calculate the molecular weight of ‘CH4’
We know the atomic weight of ‘C’ and ‘H’.
We will use those values to calculate the molecular weight.
- Getting the no. of characters inside that string and then comparing each character with all the basic elements and if it is found than assigning the value of atomic weight to it.
- It means each character will be searched according to no. of elements in there.
- Now is the case if the character is a no. itself than we have to add that character atm weight value the no. of times the no. which has been there inside the string.
- In our case as it is CH4, the above statement implies that mol. Weight of
CH4=C+H+H+H+H
- Now the length of string was 3 means it will give us the value of CH2 and 2 values will get remained, over here what we need to do is to whatever the value is inside the string we need to subtract it by 2 and multiply it with the previous character value and add it to the original weight.
- As the no. should not iterate the times of elements present, we use ‘break’ to end the inner loop once the value is being calculated.
- This is how values of each and every species will be calculated and results are shown below.
(Note: as we are searching by each element ‘AR’ cannot be searched therefore keeping it in last column so that if its weight comes out to be zero it can be assigned with the help of last value of inner loop.)
As the main progaram should give only the plots and the molecular values , commenting down the saving commands and the whole program is shown below:
close all
clear
clc
%% Opening The 'THERMO.dat' File
f1=fopen('THERMO.dat','r');
% Universal gas constant
R=8.314;
%% To Read The Header Information
for i=1:5
if i==2
tempratures=fgetl(f1);
temp=strsplit(tempratures,'');
Tmp=str2num(temp{1});
lower_global_tmp=Tmp(1);
mid_gloabal_tmp=Tmp(2);
higher_gloabal_tmp=Tmp(3);
else
comments{i}=fgetl(f1());
end
end
%% Storing the data of each species in a single cell matrix
for j=1:53
for k=1:4
A{j,k}=fgetl(f1);
end
end
%% Seprating The Sepcies coloumn and the coefficients coloumn
coefficients=A(:,2:4); % coefficeints columns
num_values=coefficients'; %transposing the cell so that each species coffiecient
%appears in a single row
species=A(:,1); % getting the species column
%converting it into char array
sp=char(species);
%% evaluating the local temperatures as well as creating dividing the rows of each species
for m=1:53
% getting the local tempratures
values{m}=strsplit(sp(m,1:end-1),'G ');
L_tmp=values{m};
local_tmp(m)=L_tmp(1,2);
%for creating species
elements_split{m}=strsplit(sp(m,:),' ');
each_species=elements_split{m};
Species(m)=each_species(1);
%seprating coefficients of each line for every species
w=num_values(:,m);
line1_species(m)=w(1,1);
line2_species(m)=w(2,1);
line3_species(m)=w(3,1);
end
L_Temp=char(local_tmp);
local_tempreatures=str2num(L_Temp);
lower_local_temp=local_tempreatures(:,1);
higher_local_temp=local_tempreatures(:,2);
mid_local_temp=local_tempreatures(:,3);
%% extracting cofficents from lines of every species
%converting it into character form so that we can use
%the function strfind
line1=char(line1_species);
line2=char(line2_species);
line3=char(line3_species);
for s=1:53
a=strfind(line1(s,:),'E');
% upper tmperature cofficients
Ah1(s,:)=line1(s,1:a(1)+3);
ah1(s)=str2double(Ah1(s,:));
Ah2(s,:)=line1(s,a(1)+4:a(2)+3);
ah2(s)=str2double(Ah2(s,:));
Ah3(s,:)=line1(s,a(2)+4:a(3)+3);
ah3(s)=str2double(Ah3(s,:));
Ah4(s,:)=line1(s,a(3)+4:a(4)+3);
ah4(s)=str2double(Ah4(s,:));
Ah5(s,:)=line1(s,a(4)+4:a(5)+3);
ah5(s)=str2double(Ah5(s,:));
Ah6(s,:)=line2(s,1:a(1)+3);
ah6(s)=str2double(Ah6(s,:));
Ah7(s,:)=line2(s,a(1)+4:a(2)+3);
ah7(s)=str2double(Ah7(s,:));
%lower temperature coffiecients
AL1(s,:)=line2(s,a(2)+4:a(3)+3);
aL1(s)=str2double(AL1(s,:));
AL2(s,:)=line2(s,a(3)+4:a(4)+3);
aL2(s)=str2double(AL2(s,:));
AL3(s,:)=line2(s,a(4)+4:a(5)+3);
aL3(s)=str2double(AL3(s,:));
AL4(s,:)=line3(s,1:a(1)+3);
aL4(s)=str2double(AL4(s,:));
AL5(s,:)=line3(s,a(1)+4:a(2)+3);
aL5(s)=str2double(AL5(s,:));
AL6(s,:)=line3(s,a(2)+4:a(3)+3);
aL6(s)=str2double(AL6(s,:));
AL7(s,:)=line3(s,a(3)+4:a(4)+3);
aL7(s)=str2double(AL7(s,:));
% aligning the contents of each species
c(s,:)=[ah1(s) ah2(s) ah3(s) ah4(s) ah5(s) ah6(s) ah7(s) aL1(s) aL2(s) aL3(s) aL4(s) aL5(s) aL6(s) aL7(s)];
end
%% calucaltion for Cp,H,S
% Cp/R=a1+(a2*T)+(a3*T^2)+(a4*T^3)+(a5*T^4)
% H/R=(a1*T)+((a2*T^2)/2)+((a3*T^3)/3)+((a4*T^4)/4)+((a5*T^5)/5)+a6
% S/R=a1*lnT + a2*T + ((a3*T^2)/2) + ((a4*T^3)/3)+((a5*T^4)/4)+a7
for p=1:53
T=lower_local_temp(p):higher_local_temp(p);
if T<=mid_local_temp(p)
Cp=R*(aL1(p)+(aL2(p)*T)+(aL3(p)*(T.^2))+(aL4(p)*(T.^3))+(aL5(p)*(T.^4)));
H=R*((aL1(p)*T)+((aL2(p)*(T.^2))/2)+((aL3(p)*(T.^3))/3)+((aL4(p)*(T.^4))/4)+((aL5(p)*(T.^5))/5)+aL6(p));
S=R*((aL1(p)*log(T)) + (aL2(p)*T) + ((aL3(p)*(T.^2))/2) + ((aL4(p)*(T.^3))/3)+((aL5(p)*(T.^4))/4)+(aL7(p)));
else
Cp=R*(ah1(p)+(ah2(p)*T)+(ah3(p)*(T.^2))+(ah4(p)*(T.^3))+(ah5(p)*(T.^4)));
H=R*((ah1(p)*T)+((ah2(p)*(T.^2))/2)+((ah3(p)*(T.^3))/3)+((ah4(p)*(T.^4))/4)+((ah5(p)*(T.^5))/5)+ah6(p));
S=R*((ah1(p)*log(T)) + (ah2(p)*T) + ((ah3(p)*(T.^2))/2) + ((ah4(p)*(T.^3))/3)+((ah5(p)*(T.^4))/4)+(ah7(p)));
end
figure(p)
plot(T,Cp,'r-','linewidth',2);
title(strcat('Specific Heat Vs Temperature of:',Species{:,p}))
xlabel('Temperature(K)')
ylabel('specific Heat(J/mol.K)')
grid on
% foldername=sprintf('%s',Species{:,p});
%
% mkdir(foldername(1,:))
% picturename = fullfile(foldername(1,:), 'Specific_Heat');
% saveas (figure(p), picturename, 'png')
figure(p+53)
plot(T,H,'b-','linewidth',2);
title(strcat('Enthalpy Vs Temperature of:',Species{:,p}))
xlabel('Temperature(K)')
ylabel('Enthalpy(J)')
grid on
% picturename = fullfile(foldername(1,:), 'Enthaply');
% saveas (figure(p+53), picturename, 'png')
figure(p+106)
plot(T,S,'g-','linewidth',2);
title(strcat('Entropy Vs Temperature of:',Species{:,p}))
xlabel('Temperature(K)')
ylabel('Entropy(J/K)')
grid on
% picturename = fullfile(foldername(1,:), 'Entropy');
% saveas (figure(p+106), picturename, 'png')
end
%% calculation of molecular weight
aw=[16.00,1.00,12.00,14.00,40.00];
elements= ["O","H","C","N","AR"];
for i=1:53
z=char(Species(i));
weight=0;
for j=1:length(z)
for k=1:length(elements)
if strcmpi(z(j),elements(k))
weight=weight+aw(k);
r=k;
else
if str2num(z(j))>=1
val=str2num(z(j));
weight=weight+aw(r)+(aw(r)*(val-2));
break;
end
end
end
end
%specially for AR
if weight==0
weight=aw(k);
end
wt(i)=weight;
disp(strcat('Molecular weight of:',Species{i},'=',num2str(wt(i))));
end
PLOTS which were needed are shown below and molecular weight is already been shown above.
1)O2
2)N2
3)CO2
CONCLUSION:
The graphs of specific heat , entropy and enthalpy has been plotted with the use of ‘THERMO.dat’ file and molecular weight has also being calculated.