Challenge

Project 1 - Parsing NASA thermodynamic data

AAnu Anilkumar2 Feb 2021

Objective :

1) To write a function to extract 14 co-efficients from the given data file and calculate Specific heat (Cp) , Enthalpy (H)

and entropy(S) for all the species .

2) To calculate Molecular Weight for all species and display it .

3) Plot Cp , H and S vs Temperature for local temperature range for all species .

4) Save these plots with suitable names and dump them in seperate folders .

5)To write detailed report for the codes .

Theory :

a) Specific Heat Capacity (Cp)

Specific heat is the quantity of heat required to raise temperature of 1kg of a substance by 1 degree Celcius.

b)Enthalpy (H)

Enthalpy is defined as the sum of internal energy and product of pressure and volume of a thermodynamic system .

c)Entropy (S)

Entropy is the measure of amount of a systems thermal energy to unit temperature that is unavailable to do useful work . Since the amount of work done is obtained from ordered molecular motion , Entropy is also a measure of molecular disorder or randomness of a system .

The given task is obtain value for Cp ,H and S using the below equation .

Where R is the Universal gas constant .

Program Logic :

f1=fopen('THERMO.dat','r') : To open file THERMO.dat in read format .

for i=1:5 : For loop to skip first 5 lines of the file

fgetl(f1); : To read the particular line of file excluding newline characters

end : End of for loop

for j=1:53 : Initialising for loop to read all specie data set .

Cp=0; : Definition of Cp to change value to 0 at end of every loop
H=0; :Definition of H to change value to 0 at end of every loop
S=0; :Definition of S to change value to 0 at end of every loop
tline = fgetl(f1) : To read first specie line from the file
A1=strsplit(tline,' ') : To split the string depending on the number of spaces
specie =A1(1) : To designate specie as first cell in cell array A1 .
specie_name=char(specie) : To change string (specie) to character data type
temp_low= str2num(A1{end-3}) : To define local temperature low as 3rd cell from end of cell array A1
temp_mid=str2num(A1{end-1}) :To define local temperature mid as 1st cell from end of cell array A1
temp_high=str2num(A1{end-2}) :To define local temperature high as 2nd cell from end of cell array A1
coeff_specie1=fscanf(f1,'%f') : To obtain all values of co-efficeints as a column array of floating point numbers format and to bring the file pointer to beginning of next specie line

R=8.31446261815324 : Defining R value
Temperature = (temp_low:temp_high); : Defining Temperature array varying from local low Temp to local high temp

% To calculate Cp, H and S
a1=coeff_specie1(1); : To designate a1 as 1st element of coefficient matrix .
a2=coeff_specie1(2); :To designate a2 as 2nd element of coefficient matrix .
a3=coeff_specie1(3); :To designate a3 as 3rd element of coefficient matrix .
a4=coeff_specie1(4); :To designate a4 as 4th element of coefficient matrix . a5=coeff_specie1(5); :To designate a5 as 5th element of coefficient matrix .
a6=coeff_specie1(7); :To designate a6 as 6th element of coefficient matrix .
a7=coeff_specie1(8); :To designate a7 as 7th element of coefficient matrix .
a8=coeff_specie1(9); :To designate a8 as 8th element of coefficient matrix .
a9=coeff_specie1(10); :To designate a9 as 9th element of coefficient matrix .
a10=coeff_specie1(11); :To designate a10 as 10th element of coefficient matrix .
a11=coeff_specie1(13); :To designate a11 as 11th element of coefficient matrix .
a12=coeff_specie1(14); :To designate a12 as 12th element of coefficient matrix .
a13=coeff_specie1(15); :To designate a13 as 13th element of coefficient matrix .
a14=coeff_specie1(16); :To designate a14 as 14th element of coefficient matrix .
a=[a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14]; :To designate [a] as an array of coefficients .

k = 1; : Defining k for varying Cp,H and S

for i=(temp_low:temp_high) : For loop to run for all temperature from local low to local high

if i<temp_mid : If condition to change co-efficient value after mid temperature value
x=8; : defining x as 8 to get co-efficient value from mid to higher temp value
else : Else condition for if defined earlier
x=1; : defining x as 1 for getting co-efficient values from low to mid temp value
end : end of if-else

Cp(k)= (a(x)+a(x+1)*i+a(x+2)*i^2+a(x+3)*i^3+a(x+4)*i^4)*R; : Equation to calcuate Cp and vary co-efficient depending upon temp range

H(k)= (a(x)+a(x+1)*i/2+a(x+2)*i^2/3+a(x+3)*i^3/4+a(x+4)*i^4/5+a(x+5)/i)*R*i; : Equation to calcuate H and vary co-efficient depending upon temp range

S(k)= (a(x)*log(i) +a(x+1)*i + a(x+2)*i^2/2 + a(x+3)*i^3/3 + a(x+4)*i^4/4 + a(x+6))*R ; :Equation to calcuate S and vary co-efficient depending upon temp range

k=k+1; : Incrementing value of k to obtain new co-efficients

end : End of for loop to calculate Cp,H and S

% To calculate molecular wieght
elements = ["H","C","N","O","S","A"] ; : Defining an array of elements with element name
atomic_mass=[1.00797,12.011,14.0067,15.9994,32.06,39.948]; : defining an array of atomic mass in order of elements in array defined above
molecular_weight=0; : defining a variable molecular weight

for i= 1:length(specie_name) : For loop to check every character/number in character array specie_name

[num,status] = str2num(specie_name(i)); : using str2num function to check for numerical suffixes in specie name

if(status) : if condition to calculate molecular mass
molecular_weight=molecular_weight+(atomic_mass(n)*(num-1)); : equation to calculate molecular mass if status=1 ,
else
for k=1:length(elements) : for loop to calculate molecular mass according to else condition
if specie_name(i)== elements(k) : if condition to check if elements in specie_name is equal to elements in element matrix
molecular_weight=molecular_weight+atomic_mass(k) : equation to calculate molecular weight if status=0

n=k; : defining n to calculate atomic mass of elements with suffixes.
end :end of 2nd if condition
end :end of 2nd for loop (k)
end : end of first if condition
end :end of first for loop(i)
MOLECULAR_WEIGHT= molecular_weight ; : defining MOLECULAR_WEIGHT for diplaying the value .
fprintf('molecular weight of %s : %f \n', specie_name,MOLECULAR_WEIGHT) : To display file according to element name and molecular mass

% To plot graphs for Specific Heat , Enthalpy , Entropy
mkdir (specie_name) : command to create folder according to specie name

figure(1)
plot(Temperature,Cp,'linewidth',3) : to plot Cp Vs temperature with linewidth 3
title(' Specific Heat Vs Temperature') : To assign title to the plot "Specific heat VS Temperature"
xlabel('Temperature T') : to label x axis as Temperature T
ylabel('Specific heat Cp ') : to label y axis as Specific heat Cp
url1 = sprintf('E:/anu_matlab/project1/%s/ Specific Heat Plot.jpg', specie_name) : to assing file path and name to file according to specie name
saveas(gcf,url1) :to save the plot in url1
figure(2)
plot(Temperature,H,'linewidth',3,'color','y') :to plot Cp Vs temperature with linewidth 3 and color yellow
title(' Enthalpy Vs Temperature') :To assign title to the plot "Enthalpy Vs temperature"
xlabel('Temperature T') :to label x axis as Temperature T
ylabel('Enthalpy H') :to label y axis as Enthalpy H

url2 = sprintf('E:/anu_matlab/project1/%s/Enthalpy VS Temperature.jpg', specie_name):to assing file path and name to file according to specie name

saveas(gcf,url2) :to save the plot in url2
figure(3)
plot(Temperature,S,'linewidth',3,'color','g') :to plot Cp Vs temperature with linewidth 3 and color green
title(' Entropy Vs Temperature') :To assign title to the plot "Entropy Vs temperature"
xlabel('Temperature T') :to label x axis as Temperature T
ylabel('Entropy S') :to label y axis as Entropy S
url3 = sprintf('E:/anu_matlab/project1/%s/Entropy VS Temperature.jpg', specie_name) :to assing file path and name to file according to specie name
saveas(gcf,url3) :to save the plot in url3

end : End of for loop (j)

Programe code :

clear all 
close all
clc

% To read file THERMO.dat
f1=fopen('THERMO.dat','r')



for i=1:5
   
fgetl(f1); 
    
end
for j=1:53
    
    Cp=0;
    H=0;
    S=0;
tline = fgetl(f1)
A1=strsplit(tline,' ')
specie =A1(1)
specie_name=char(specie)
temp_low= str2num(A1{end-3})
temp_mid=str2num(A1{end-1})
temp_high=str2num(A1{end-2})
 coeff_specie1=fscanf(f1,'%f')
 R=8.31446261815324
 Temperature = (temp_low:temp_high);
 
 % To calculate Cp, H and S
 a1=coeff_specie1(1);
 a2=coeff_specie1(2);
 a3=coeff_specie1(3);
 a4=coeff_specie1(4);
 a5=coeff_specie1(5);
 a6=coeff_specie1(7);
 a7=coeff_specie1(8);
 a8=coeff_specie1(9);
 a9=coeff_specie1(10);
 a10=coeff_specie1(11);
 a11=coeff_specie1(13);
 a12=coeff_specie1(14);
 a13=coeff_specie1(15);
 a14=coeff_specie1(16);
a=[a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14];
 k = 1;
 
 for i=(temp_low:temp_high)
     
     if i<temp_mid 
         x=8;
     else
         x=1;
     end
         
 Cp(k)= (a(x)+a(x+1)*i+a(x+2)*i^2+a(x+3)*i^3+a(x+4)*i^4)*R;
 
 H(k)=  (a(x)+a(x+1)*i/2+a(x+2)*i^2/3+a(x+3)*i^3/4+a(x+4)*i^4/5+a(x+5)/i)*R*i;
 
 S(k)= (a(x)*log(i) +a(x+1)*i + a(x+2)*i^2/2 + a(x+3)*i^3/3 + a(x+4)*i^4/4 + a(x+6))*R ; 
 
 k=k+1;
 end
 
 % To calculate molecular wieght
 elements = ["H","C","N","O","S","A"] ;
 atomic_mass=[1.00797,12.011,14.0067,15.9994,32.06,39.948];
 molecular_weight=0;

 for i= 1:length(specie_name)
 
     [num,status] = str2num(specie_name(i));
     if(status)
         molecular_weight=molecular_weight+(atomic_mass(n)*(num-1));
     else   
         for k=1:length(elements)
   if  specie_name(i)== elements(k)
       molecular_weight=molecular_weight+atomic_mass(k)
       n=k;
   end
         end
     end
 end
 MOLECULAR_WEIGHT= molecular_weight ;
 fprintf('molecular weight of %s : %f \n', specie_name,MOLECULAR_WEIGHT)      
       
 
%  To create folder according to specie name
 mkdir (specie_name)
 %  To plot graphs for Specific Heat , Enthalpy , Entropy 
 figure(1)
 plot(Temperature,Cp,'linewidth',3)
 title(' Specific Heat Vs Temperature')
 xlabel('Temperature T')
 ylabel('Specific heat Cp ')
 url1 = sprintf('E:/anu_matlab/project1/%s/ Specific Heat Plot.jpg', specie_name)
 saveas(gcf,url1)
 figure(2)
 plot(Temperature,H,'linewidth',3,'color','y')
 title(' Enthalpy Vs Temperature')
 xlabel('Temperature T')
 ylabel('Enthalpy H')
 url2 = sprintf('E:/anu_matlab/project1/%s/Enthalpy VS Temperature.jpg', specie_name)
 saveas(gcf,url2)
 figure(3)
 plot(Temperature,S,'linewidth',3,'color','g')
 title(' Entropy Vs Temperature')
 xlabel('Temperature T')
 ylabel('Entropy S')
 url3 = sprintf('E:/anu_matlab/project1/%s/Entropy VS Temperature.jpg', specie_name)
 saveas(gcf,url3)
 

end

Results :

1) The 14 co-efficients for each species were obtained from the file and value of Cp , H and S were calculated.

2) Molecular weight of each specie were calculated and displayed in the command window .