Project 1 - Parsing NASA thermodynamic data

Aim:1.To write a function to extract 14 co-efficients and program for calculation of specific heat,entropy,enthalpy from given file.

2.Calculate molecular weight of each species.

3.Plot a graph for specific heat,enthalpy,entropy for temperature for each species.

Theory:File parsing in matlab in general ,the term parsing refers to MATLAB reading an m-file(script or function) and translate it into an internal code that it stores in memory.It is 'parsed' file in memory that is actually run when execute the code.

Here we extract the data from thermodynamic data file and calculate the enthalpy,entropy specific heat for all the species in file and also calculate molecular weight.

1.Enthalpy:It is a thermodynamic property,and it is the sum of internal energy of system and product of its pressure and volume.It is denoted by "H".

H=E+PV

where, E=internal energy

P=pressure

V=volume

H=enthalpy

but in terms of NASA THERMODYNAMIC POLYNOMIAL FORMULATE AS:

H/RT = a1 + a2 T /2 + a3 T^2 /3 + a4 T^3 /4 + a5 T^4 /5 + a6/T

2.Entropy:It is a thermodynamic function,it is the measure of system thermal energy per unit temperature that is unavailable for doing useful work.It is denoted by "S"

S=q/T

where S=entropy

q=thermal energy

T=temperature.

but in terms of NASA THERMODYNAMIC POLYNOMIAL FORMULATE AS:

S/R = a1 lnT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7

3.Specific heat:It is a thermodynamic property,the quantity of heat required to raise the temperature of one gram of substance by one celcius degree. It is denote by C_p.

C_p = Q/m*▲T

where, Q=heat energy

m=mass

▲T=change in temperature

but in terms of NASA THERMODYNAMIC POLYNOMIAL FORMULATE AS:

Cp/R = a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4

4.Molecular weight:The weight in atomic mass units of each of the atoms in a given formula is the molecular weight of a substance.

List of 53 Elements contained in 'THERMO.dat' file:
AR	C	C2H	C2H2	C2H3	C2H4
C2H5	C2H6	C3H7	C3H8	CH	CH2
CH2(S)	CH2CHO	CH2CO	CH2O	CH2OH	CH3
CH3CHO	CH3O	CH3OH	CH4	CN	CO
CO2	H	H2	H2CN	H2O	H2O2
HCCO	CCOH	HCN	HCNN	HCNO	HCO
HNCO	HNO	HO2	HOCN	N	N2
N2O	NCO	NH	NH2	NH3	NNH
NO	NO2	O	O2	OH

code of all functions:

1.specific heat

Cp/R = a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4

% Specific heat function
function Cp = Specificheat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp) 
% Comparing with global temperature
% If less than global temperature - minimum coefficients are taken
if Temp <= global_medium_temp
Cp = R*(a8 + (a9*Temp) + (a10*(Temp).^2) + (a11*(Temp).^3) + (a12*(Temp).^4)); 
else
% Else maximum coefficients are taken
Cp = R*(a1 + (a2*Temp) + (a3*(Temp).^2) + (a4*(Temp).^3) + (a5*(Temp).^4));
end

2.Entropy

S/R = a1 lnT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7

% Entropy function
function S = Entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,Temp,R,global_medium_temp) 
% Comparing with global temperature
% If less than global temperature - minimum coefficients are taken
 if Temp <= global_medium_temp
 S = R.*(a8.*log(Temp) + (a9.*Temp) + (a10*((Temp).^2)/2) + (a11*((Temp).^3)/3) + (a12*((Temp).^4)/4) + a14);
 else
% Else - maximum coefficients are taken
 S = R.*(a1.*log(Temp) + (a2.*Temp) + (a3*((Temp).^2)/2) + (a4*((Temp).^3)/3) + (a5*((Temp).^4)/4) + a7);
end

3.Enthalpy

H/RT = a1 + a2 T /2 + a3 T^2 /3 + a4 T^3 /4 + a5 T^4 /5 + a6/T

% Enthalpy function
function H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,Temp,R,global_medium_temp) 
% Comparing with global temperature
% If less than global temperature - minimum coefficients are taken
 if Temp <= global_medium_temp
 H = R*Temp.*(a8 + ((a9*Temp)/2) + ((a10*(Temp).^2)/3) + ((a11*(Temp).^3)/4) + ((a12*(Temp).^4)/5) +((a13./Temp)));
 else
% Else maximum coefficients are taken 
 
 H = R*Temp.*(a1 + ((a2*Temp)/2) + ((a3*(Temp).^2)/3) + ((a4*(Temp).^3)/4) + ((a5*(Temp).^4)/5) +((a6./Temp)));
end

4.Molecular weight

% Molecular Weight function
function Molecularweight = Molecular_weight(Species_SP)
 % Standard - Atomic Name
 atomic_name = ['H','C','O','N','A']; 
 % Hydrogen, Carbon, Oxygen, Nitrogen, Argon
 
 fprintf('Molecular Weight of ');
 fprintf(Species_SP);
 fprintf(' : ');
 % Standard Atomic Number
 Atomic_Weight = [1.008 12.011 15.999 14.007 39.948]; 
 % Starting Range of Molecule 
 Molecularweight = 0; 
 
 % Adding Atomic number to obtain particular value
 for i = 1:length(Species_SP)
 for j = 1:length(atomic_name)
 if strcmpi(Species_SP(i),atomic_name(j))
 Molecularweight = Molecularweight +Atomic_Weight(j);
 position = j;
 end
 end
 n = str2num(Species_SP(i));
 
 if n>1
 
 Molecularweight = Molecularweight + Atomic_Weight(position)*(n-1);
 
 end
 end
 fprintf('molecular weight %s :',Species_SP);
 fprintf('%d',Molecularweight);
 disp(' '); 
end

code explanation:

first of all we will make path for program to extract the all data of thermo.dat file.

after that first confirm that your thermo.dat file save in same path of program folder

overall code:

clear all
close all
clc
% Make directory
mkdir 'D:\pratik.tembhare(skill lync)'
% Reading Data file - NASA Thermodynamic
% Open a file & 'r' Open the file for reading 
f1=fopen('THERMO.dat','r')
gt=fgetl(f1);
% Obtaing Temperatures - by Reading Header Information
global_temp=fgetl(f1);
% Spliting
t=strsplit(global_temp);
% Converting string cells to numbers using str2num,
global_low_temp=str2num(t{2});
global_medium_temp=str2num(t{3});
global_high_temp=str2num(t{4});
% Temperature range & Gas Constant
 % Temperature Range assumed
Temp=linspace(global_low_temp,global_high_temp,1000);
% Universal Gas Constant J/(mol-K)
 R = 8.314; 
 
% skipping three comment lines in the file
for i = 1:3
 gt=fgetl(f1);
end
 for i = 1:53
gt=fgetl(f1);
A=strsplit(gt,' ');
Species_SP=(A{1})
line_1=fgetl(f1) ; 
a=strfind(line_1,'E');
% Extracting higher and lower temperature coefficients
% % Extracting FIRST 7 coefficients i.e HIGH-temperature coefficients 
a1=line_1(1:a(1)+3);
a1=str2num(a1); %1
a2=line_1(a(1)+4:a(2)+3);
a2=str2num(a2) ; %2
a3=line_1(a(2)+4:a(3)+3);
a3=str2num(a3) ; %3
a4=line_1(a(3)+4:a(4)+3);
a4=str2num(a4); %4
a5=line_1(a(4)+4:a(5)+3);
a5=str2num(a5) ; %5
line_2=fgetl(f1);
b=strfind(line_1,'E');;
a6=line_2(1:b(1)+3);
a6=str2num(a6) ; %6
a7=line_2(a(1)+4:b(2)+3);
a7=str2num(a7); %7
 % Extracting SECOND 7 coefficients - LOW-temperature coefficients
a8=line_2(b(2)+4:b(3)+3); 
a8=str2num(a8); %8
a9=line_2(b(3)+4:b(4)+3);
a9=str2num(a9) ; %9
a10=line_2(b(4)+4:b(5)+3);
a10=str2num(a10) ; %10
line_3=fgetl(f1);
c=strfind(line_1,'E');
a11=line_3(1:c(1)+3);
a11=str2num(a11); %11
a12=line_3(c(1)+4:c(2)+3);
a12=str2num(a12); %12
a13=line_3(c(2)+4:c(3)+3);
a13=str2num(a13); %13
a14=line_3(c(3)+4:c(4)+3);
a14=str2num(a14); %14
% Specific heat calculation
 Cp = Specificheat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,Temp,R,global_medium_temp);
 % Enthalpy calculation
 H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,Temp,R,global_medium_temp);
 % Entropy calculation
 S = Entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,Temp,R,global_medium_temp);
 % Calculating molecular weight of species
 Molecularweight(i) = Molecular_weight(Species_SP);
 
 % Read & Write Molecular list
 Molecular_list = fopen('Molecular weight.txt','a+');
 fprintf( Molecular_list,'Molecular weight of %s is %d n',Species_SP,Molecularweight);
 fclose(Molecular_list);
 % plotting(CP,H,S,Temp,Species_SP)
 % Saving as a folder
 cur_dir=pwd;
 mkdir(Species_SP);
 cd(Species_SP);
 % Plot 1, Specific heat (CP) vs Temperature Range (Temp)
 figure(1)
 plot(Temp,Cp,'linewidth',2,'color','r');
 xlabel('Temperature [K]');
 ylabel('Specific Heat [kJ]');
 title(sprintf('Varying Specific Heat with Temperature Range of %s',Species_SP));
 grid on
 saveas(figure(1),'Specific Heat.jpg');
 % Plot 2, Enthalpy (H) vs Temperature Range (Temp)
 figure(2)
 plot(Temp,H,'linewidth',2,'color','g');
 xlabel('Temperature [K]');
 ylabel('Enthalpy in [kJ/(mol)]');
 title(sprintf('Varying Enthalpy with Temperature Range of %s',Species_SP));
 grid on
 saveas(figure(2),'Ethalpy.jpg');
 % Plot 3, Entropy (S) vs Temperature Range (Temp)
 figure(3)
 plot(Temp,S,'linewidth',2,'color','b');
 xlabel('Temperature [K]');
 ylabel('Entropy in [kJ/(mol-K)]');
 title(sprintf('Varying Entropy with Temperature Range of %s',Species_SP));
 grid on
 saveas(figure(3),'Entropy.jpg');
 cd(cur_dir);
 %MOLECULAR WEIGHT
 M = Molecular_weight(Species_SP)
 f2 = fopen('molecular mass','w');
 fprintf('%s n %s n %f n','species','Molecular mass',M)
 fclose(f2)
end

output: