To obtain the plots of enthalpy, entropy and specific heat with temperature of various compounds by parsing the given thermodynamic data using Matlab/Octave programming
AIM
To obtain the plots of enthalpy, entropy and specific heat with temperature of various compounds by parsing the given thermodynamic data using Matlab/Octave programming.
PARSING
Parsing data means breaking it down into components or parts so that its syntax and parameters can be divided and understood better. The parsing tool takes the tool in the form of a file or program instructions and builds a data structure in the form of a abstract syntax tree. The overall process of parsing involves three stages.
1) lexical analyis - Tokens are produced from a stream of input characters which are then broken into smaller components.
2) syntactic analysis - Checks whether the tokens generated form a meaningful expression or not.
3) semantic analysis - The meaning of the validated expression are determined .
Parsers are widely used in programming languages such as java, html,C, C++ and other modelling and scripting languages.
GOVERNING EQUATIONS
`(Cp)/R=a1+a2*T+a3*T^2+a4*T^3+a5*T^4`
`H/(RT)=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,
Cp = specific heat
H = enthalpy
S = entropy
a1 to a14 are the co-efficients
OBJECTIVE
The main objective of this project is to extract the 14 co-efficients as given in the thermodynamic data file provided and calculate the parameters as specific heat, enthalpy and entropy for the various compounds. The variations of these parameters with respect to temperature are plotted and the various graphs are hence obtained.
PROCEDURE
1) Open the file containing the data using the fopen command and read the lines which are only required and omit the unnecesaary ones using the fgetl command.
2) Extract the global temperature ranges using the fgetl and strsplit commands.
3) Calculate and extract all 14 co-efficients required.
4) Using for loop, run the loop for all plots and graphs for all compounds after calculating the values of all the co-efficients.
5) Calculate the specific heat, enthalpy and entropy
6) Store all plots of a compound in a file and store all such files.
7) Calculate the molecular weights.
CODE EXPLANATION
The code consists of various portions. The main code along with the functions to call the specific heat, enthalpy and entropy are also mentioned.
1) Main code
f1 = fopen('THERMO.dat','r'); - This statement is used to open the file ‘THERMO.DAT’ for studying its data.
header = fgetl(f1); - This statement reads the first line of the file which is actually of no use in the program.
global_temp = fgetl(f1); - This statement reads the next line where there are the ranges of global temperatures.
S = strsplit(global_temp,' '); - The ‘strsplit’ command is used to identify the various elements in a string and separate them out wherever there is a space in between two successive values/words.
global_temp_low = str2num(S{2}); - Used to extract the lower global temperature.
global_temp_mid = str2num(S{3}); - Used to extract the mid global temperature.
global_temp_high = str2num(S{4});- Used to extract the higher global temperature.
fgetl(f1); - Used to read the next line which is ignored.
fgetl(f1); - Used to read the next line which is ignored.
fgetl(f1); - Used to read the next line which is ignored.
for i = 1:53 – To read the various compounds values from 1 to 53.
L = fgetl(f1); - Used to read line by line of the values.
A = strsplit(L,' ' ); - It is used to separate the values as individual ones wherever there is a space between two values.
spec = A(1); - Each individual value is stored in a separate variable.
M = strfind(L,'G'); - The obtainment of values is based on the location of the ‘G’ in them. In order to call the respective values later, the strfind option finds where all the letter is present and gives the location.
local_temp_low = L(M+3:55); - Gets the lower local temperature value.
local_temp_mid = L(M+13:65); - Gets the mid local temperature value
local_temp_high= L(M+23:75); - Gets the higher local temperature value
ltl = str2double(local_temp_low); - Stores the value in a separate variable as a double/numeric value so that calculations can be performed.
ltm = str2double(local_temp_mid); Stores the value in a separate variable as a double/numeric value so that calculations can be performed.
lth = str2double(local_temp_high); Stores the value in a separate variable as a double/numeric value so that calculations can be performed.
line1 = fgetl(f1);
X = strfind(line1,'E');
a1 = line1(1:X(1)+3);
a1 = str2double(a1);
a2 = line1(X(1)+4:X(2)+3);
a2 = str2double(a2);
a3 = line1(X(2)+4:X(3)+3);
a3 = str2double(a3);
a4 = line1(X(3)+4:X(4)+3);
a4 = str2double(a4);
a5 = line1(X(4)+4:X(5)+3);
a5 = str2double(a5); - The following lines are used to extract the first five co-efficients based on the location of the letter ‘E’ in the particular read line.
line2 = fgetl(f1);
Y = strfind(line2,'E');
a6 = line2(1:Y(1)+3);
a6 = str2double(a6);
a7 = line2(Y(1)+4:Y(2)+3);
a7 = str2double(a7);
a8 = line2(Y(2)+4:Y(3)+3);
a8 = str2double(a8);
a9 = line2(Y(3)+4:Y(4)+3);
a9 = str2double(a9);
a10 = line2(Y(4)+4:Y(5)+3);
a10 = str2double(a10); - The following lines are used to extract the next five co-efficients based on the location of the letter ‘E’ in the particular read line.
line3 = fgetl(f1);
Z = strfind(line3,'E');
a11 = line3(1:Z(1)+3);
a11 = str2double(a11);
a12 = line3(Z(1)+4:Z(2)+3);
a12 = str2double(a12);
a13 = line3(Z(2)+4:Z(3)+3);
a13 = str2double(a13);
a14 = line3(Z(3)+4:Z(4)+3);
a14 = str2double(a14); - The following lines are used to extract the last four co-efficients based on the location of the letter ‘E’ in the particular read line.
R = 8.314; - Initializing the constant value of gas constant.
Temp = linspace(300,5000,100); - Varying temperature from the global low limit to the global high limit with 100 equal values in between.
Cp = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Temp,R,ltl,ltm,lth); - Calling the function from the specific heat program in order to obtain the values as per the formulae.
S = entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Temp,R,ltl,ltm,lth); - Calling the function from the entropy program in order to obtain the values as per the formulae.
H = enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Temp,R,ltl,ltm,lth); - Calling the function from the enthalpy program in order to obtain the values as per the formulae.
Spec = char(spec); - Used to get the character value of the values stored.
fd = sprintf('%s',Spec); - The sprint command is used to return the output as a string.
C = pwd;
mkdir(Spec); - The mkdir command allows the creation of directories and set the permission for the directories as well.
cd(Spec); - the cd command is used to change the current working directory in operating systems. It is also used to shell scripts and batch files.
figure(1)
plot(Temp,Cp,'g');
xlabel('Temperature in Kelvin');
ylabel('Specific Heat in Kilojoules');
title('Specific heat vs Temperature'); - Plotting the graphs of temperature vs specific heat for all compounds as all this is inside the for loop.
figure(2)
plot(Temp,H,'b');
xlabel('Temperature in Kelvin');
ylabel('Enthalpy in Kilojoules/mol');
title('Enthalpy vs Temperature'); - Plotting the graphs of temperature vs enthalpy for all compounds as all this is inside the for loop.
figure(3)
plot(Temp,S,'r');
xlabel('Temperature in Kelvin');
ylabel('Entropy in Kilojoules/mol K');
title('Entropy heat vs Temperature'); - Plotting the graphs of temperature vs entropy for all compounds as all this is inside the for loop.
saveas (figure(1),'Cp.png'); - The saveas option in octave is used to save each of the plot obtained as a figure in png format with the name ‘Cp’ for all compounds. The command says to save each of the first outputs as ‘Cp’.
saveas(figure(2),'H.png'); - The saveas option in octave is used to save each of the plot obtained as a figure in png format with the name ‘H’ for all compounds. The command says to save each of the second outputs as ‘H’.
saveas(figure(3),'S.png'); - The saveas option in octave is used to save each of the plot obtained as a figure in png format with the name ‘S’ for all compounds. The command says to save each of the first outputs as ‘S’.
cd(C);
Mol_mass = molecular_mass(spec) – Calling the molecular masses function
f2 = fopen('Molecular Masses','w'); - Used to open the file
fprintf('%s n %s n %f n','spec','Molecular Mass','Mol_mass') – prints the molecular masses.
fclose(f2)
endfor
fclose(f1); - the fclose command is used to closed the opened file.
2) Specific heat function code
function Cp = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Temp,R,ltl,ltm,lth); - The function in octave is used to call the values in another program where typing all this there becomes a tedious process.
if(Temp>=ltm)
Cp = (a1+a2.*(Temp)+a3.*(Temp.^2)+a4.*(Temp.^3)+a5.*(Temp.^4)).*R
Else
Cp = (a8.*Temp+a9.*(Temp.^2)./2+a10.*(Temp.^3)./3+a11.*(Temp.^4)./4+a12.*(Temp.^5)./5+a13).*R
end
end
3) Enthalpy function code
function H = enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Temp,R,ltl,ltm,lth); - The function in octave is used to call the values in another program where typing all this there becomes a tedious process.
if(Temp>=ltm)
H = (a1.*Temp+a2.*(Temp.^2)./2+a3.*(Temp.^3)./3+a4.*(Temp.^4)./4+a5.*(Temp.^5)./5+a6).*R
else
H = (a8.*Temp+a9.*(Temp.^2)./2+a10.*(Temp.^3)./3+a11.*(Temp.^4)./4+a12.*(Temp.^5)./5+a13).*R
end
end
4) Entropy function code
function S = entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Temp,R,ltl,ltm,lth); - The function in octave is used to call the values in another program where typing all this there becomes a tedious process
if(Temp>=ltm)
S = ((a1.*(log(Temp))+a2.*Temp+a3.*(Temp.^2)./3+a4.*(Temp.^3)./3+a5.*(Temp.^4)./4+a7)).*R
else
S = ((a8*(log(Temp))+a9.*Temp+a10.*(Temp.^2)./2+a11.*(Temp.^3)./3+a12.*(Temp.^4./4+a14))).*R
end
5) Molecular masses function code
function Mol_mass = molecular_mass(spec) - The function in octave is used to call the values in another program where typing all this there becomes a tedious process
atom = { 'H','C','N','O','AR'}; - Initializing the values.
atom_weight = { 1 12 14 16 40 }; - Initializing the atomic weights.
Mol_mass = 0;
for i=1:length(spec)
for j=1:length(atom_weight)
if(strcmp(spec(i),atom(j))==1)
Mol_mass=Mol_mass+atom_weight{j};
g=j;
end
end
N=str2double(spec(i));
if(N>1)
Mol_mass=Mol_mass+(atom_weight(g)*(N-1));
end
PLOTS
1) FOR OXYGEN(02)
2) FOR NITROGEN(N2)
3) FOR CARBON DIOXIDE(CO2)
THE BELOW SCREENSHOT SHOWS THE VARIOUS FOLDERS SAVED AS PER THE DATA
CODE
The code consists of the main code and the functions for calling various parameters.
1) Main code
f1 = fopen('THERMO.dat','r'); % Opening the data file
header = fgetl(f1); % Reading the header
% Getting global temperature values
global_temp = fgetl(f1);
S = strsplit(global_temp,' ');
global_temp_low = str2num(S{2});
global_temp_mid = str2num(S{3});
global_temp_high = str2num(S{4});
fgetl(f1);
fgetl(f1);
fgetl(f1);
for i = 1:53
L = fgetl(f1);
A = strsplit(L,' ' );
spec = A(1);
M = strfind(L,'G');
local_temp_low = L(M+3:55);
local_temp_mid = L(M+13:65);
local_temp_high= L(M+23:75);
ltl = str2double(local_temp_low);
ltm = str2double(local_temp_mid);
lth = str2double(local_temp_high);
% Extracting first five co-efficients
line1 = fgetl(f1);
X = strfind(line1,'E');
a1 = line1(1:X(1)+3);
a1 = str2double(a1);
a2 = line1(X(1)+4:X(2)+3);
a2 = str2double(a2);
a3 = line1(X(2)+4:X(3)+3);
a3 = str2double(a3);
a4 = line1(X(3)+4:X(4)+3);
a4 = str2double(a4);
a5 = line1(X(4)+4:X(5)+3);
a5 = str2double(a5);
% Extracting next five co-efficients
line2 = fgetl(f1);
Y = strfind(line2,'E');
a6 = line2(1:Y(1)+3);
a6 = str2double(a6);
a7 = line2(Y(1)+4:Y(2)+3);
a7 = str2double(a7);
a8 = line2(Y(2)+4:Y(3)+3);
a8 = str2double(a8);
a9 = line2(Y(3)+4:Y(4)+3);
a9 = str2double(a9);
a10 = line2(Y(4)+4:Y(5)+3);
a10 = str2double(a10);
% Extracting last four co-efficients
line3 = fgetl(f1);
Z = strfind(line3,'E');
a11 = line3(1:Z(1)+3);
a11 = str2double(a11);
a12 = line3(Z(1)+4:Z(2)+3);
a12 = str2double(a12);
a13 = line3(Z(2)+4:Z(3)+3);
a13 = str2double(a13);
a14 = line3(Z(3)+4:Z(4)+3);
a14 = str2double(a14);
R = 8.314;
Temp = linspace(300,5000,100);
Cp = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Temp,R,ltl,ltm,lth);
S = entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Temp,R,ltl,ltm,lth);
H = enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Temp,R,ltl,ltm,lth);
% Changing directory and obtaining string output
Spec = char(spec);
fd = sprintf('%s',Spec);
C = pwd;
mkdir(Spec);
cd(Spec);
% Plotting temp vs Cp graphs
figure(1)
plot(Temp,Cp,'g');
xlabel('Temperature in Kelvin');
ylabel('Specific Heat in Kilojoules');
title('Specific heat vs Temperature');
% Plotting temp vs H graphs
figure(2)
plot(Temp,H,'b');
xlabel('Temperature in Kelvin');
ylabel('Enthalpy in Kilojoules/mol');
title('Enthalpy vs Temperature');
% Plotting temp vs S graphs
figure(3)
plot(Temp,S,'r');
xlabel('Temperature in Kelvin');
ylabel('Entropy in Kilojoules/mol K');
title('Entropy heat vs Temperature');
% Saving each plot
saveas (figure(1),'Cp.png');
saveas(figure(2),'H.png');
saveas(figure(3),'S.png');
cd(C);
Mol_mass = molecular_mass(spec)
f2 = fopen('Molecular Masses','w');
fprintf('%s n %s n %f n','spec','Molecular Mass','Mol_mass')
fclose(f2)
endfor
fclose(f1);2) Specific heat code
function Cp = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Temp,R,ltl,ltm,lth);
if(Temp>=ltm)
Cp = (a1+a2.*(Temp)+a3.*(Temp.^2)+a4.*(Temp.^3)+a5.*(Temp.^4)).*R
else
Cp = (a8.*Temp+a9.*(Temp.^2)./2+a10.*(Temp.^3)./3+a11.*(Temp.^4)./4+a12.*(Temp.^5)./5+a13).*R
end
end
3) Enthalpy code
function H = enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Temp,R,ltl,ltm,lth);
if(Temp>=ltm)
H = (a1.*Temp+a2.*(Temp.^2)./2+a3.*(Temp.^3)./3+a4.*(Temp.^4)./4+a5.*(Temp.^5)./5+a6).*R
else
H = (a8.*Temp+a9.*(Temp.^2)./2+a10.*(Temp.^3)./3+a11.*(Temp.^4)./4+a12.*(Temp.^5)./5+a13).*R
end
end4) Entropy code
function S = entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,Temp,R,ltl,ltm,lth);
if(Temp>=ltm)
S = ((a1.*(log(Temp))+a2.*Temp+a3.*(Temp.^2)./3+a4.*(Temp.^3)./3+a5.*(Temp.^4)./4+a7)).*R
else
S = ((a8*(log(Temp))+a9.*Temp+a10.*(Temp.^2)./2+a11.*(Temp.^3)./3+a12.*(Temp.^4./4+a14))).*R
end
5) Molecular mass code
function Mol_mass = molecular_mass(spec)
atom = { 'H','C','N','O','AR'};
atom_weight = { 1 12 14 16 40 };
Mol_mass = 0;
for i=1:length(spec)
for j=1:length(atom_weight)
if(strcmp(spec(i),atom(j))==1)
Mol_mass=Mol_mass+atom_weight{j};
g=j;
end
end
end
N=str2double(spec(i));
if(N>1)
Mol_mass=Mol_mass+(atom_weight(g)*(N-1));
end