MATLAB : Parsing thermodynamic data of NASA
AIM OF THE PROJECT
To parse NASA thermodynamic data into matlab using multiple under-mentioned objectives using the 'THERMO.dat' file as an input in which there are temperature range and 53 species, where each species contains 14 coefficients and temperature ranges. We have to parse these 14 coefficients to find Specific heat, Entropy and Enthalpy and Molecular weight.
OBJECTIVE
1. To write a function that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the data file.
2. To calculate the molecular weight of each species and display them in the command window.
3. To plot the Heat, Enthalpy and Entropy for the local temperature range (low temperature to the high temperature) specific for each species.
4. To save the plots as images with appropriate names and dump them in separate folders for each species with suitable name. All these should be generated automatically.
GOVERNING FORMULAS
1. Specific Heat :
`Cp = R*(a1 + (a2*T) + (a3*(T).^2) + (a4*(T).^3) + (a5*(T).^4))`
2. Entropy :
`S = R*(a1.*log(T) + (a2*T) + ((a3*(T).^2)/2) + ((a4*(T).^3)/3) + ((a5*(T).^4)/4) + a7)`
3. Enthalpy :
`H = R*T.*(a1 + ((a2*T)/2) + ((a3*(T).^2)/3) + ((a4*(T).^3)/4) + ((a5*(T).^4)/5) + ((a6./T)))`
where,
R = Universal Gas constant
T = Temperature (K)
Cp = Specific heat (KJ)
S = Entropy(KJ/mol-K)
H = Enthalpy(KJ/mol)
a(i) = Coefficients where, i = 1 to 14
THEORY
1. File Parsing:
Parsing, syntax analysis, or syntactic analysis is the process of analyzing a string of symbols, either in natural language, computer languages or data structures, conforming to the rules of a formal grammar. The term parsing comes from Latin pars (orationis), meaning part (of speech).Within computational linguistics the term is used to refer to the formal analysis by a computer of a sentence or other string of words into its constituents, resulting in a parse tree showing their syntactic relation to each other, which may also contain semantic and other information.Some parsing algorithms may generate a parse forest or list of parse trees for a syntactically ambiguous input.Parsing a file means reading in a data stream of some sort and building an in memory model of the semantic content of that data. This aims at facilitating perfoming some kind of transformation on the data.
2. Specific Heat:
The specific heat is the amount of heat per unit mass required to raise the temperature by one degree Celsius. The relationship between heat and temperature change is usually expressed in the form, where c is the specific heat.
3. Entropy:
Entropy, the measure of a system's thermal energy per unit temperature that is unavailable for doing useful work. Because work is obtained from ordered molecular motion, the amount of entropy is also a measure of the molecular disorder, or randomness, of a system.
4. Enthalpy:
Enthalpy, a property of a thermodynamic system, is equal to the system's internal energy plus the product of its pressure and volume. In a system contained so as to prevent mass transfer, for processes at constant pressure, the heat absorbed or released equals the change in enthalpy.
WRITING THE MAIN CODE
1. First, we will input the parameter R and will open THERMO.dat file and make this file read using fopen command.
2.Then, we will read first THERMO using fgetl command and further, read temperature file and will split the temperatures using strsplit command to convert it into cell array and next, we will convert this string to number using str2num command to find low, high and mid temperatures.
3. Next, we will assign temperature range using linspace command.
4. Now, we run loop to skip command lines the data file using fgetl command inside the loop.
5. Next, we will run the loop to find 53 element in the data file by assigning fgetl command and then spliting into the cell array to find the species line and then find the species name and then find the number of the temperatures in the species line.
6. Then, we read the new line. In order to read coefficient from the data, file we will first find the 'E' string using strfind command to find before and after numbers of the coefficients. By doing so, we will find remaining 14 coefficients and convert them to numbers from strings.
NOTE: Out of which we will find first 7 Coefficients which are for Higher temperature and then Second 7 Coefficients which are for Lower temperature.
7. Now, we will specify Specific heat Cp, Entropy S, Enthalpy H functions to find these parameters using these 14 coefficients for Lower and higher Temperature.
8.Next, we will specify molecular weight W function to find the Molecular weight of the different species.
9. Now, we will write the Molecular weight.txt file and fprintf all the molecular weights of the species.
10. And then, we will close the txt file using fclose command.
11. Finally, we change the directory using cd command and then cwe create a directory folder where we have to save our files using mkdir command and next, we will change the directory to the created folder.
12. Now, we will plot Specific heat, Entropy and Enthalpy with the Temperature for different species.
13. We we will use scanf command to in the title command to create different species title.
14. Then we will use saveas command to save the plots of different species in the directory folder.
15. Now we will change the directory to which our main code file runs in the MATLAB.
MAIN CODE:
%program for file parsing of NASA thermodynamic file
close all
clear all
clc
%opening the external file
f1= fopen('THERMO.dat','r')
%skipping this line
line_1= fgetl(f1);
%skipping this line
line_2= fgetl(f1);
%obtaining global temperature values and storing them in seperate values
global_values= strsplit(line_2,' ');
global_low_temp= str2num(global_values{2})
global_mid_temp= str2num(global_values{3})
global_high_temp= str2num(global_values{4})
%skipping this line
line_3= fgetl(f1);
%skipping the line
line_4= fgetl(f1);
%skipping the line
line_5= fgetl(f1);
%defining input universal gas constant (SI Unit: KJ/K-mol-K)
R= 8.314;
for i= 1:53
T= fgetl(f1);
%extracting the species name
species_name= T(1:9);
species_name= deblank(species_name);
%removing blank spaces
%extracting local temperatures
species_temp= T(49:75);
local_values= strsplit(species_temp,' ');
l_low_temp= str2num(local_values{1});
l_mid_temp= str2num(local_values{3});
l_high_temp= str2num(local_values{2});
%storing the local temperatures in a variable
temp= [l_low_temp l_mid_temp l_high_temp];
temp_range=linspace(l_low_temp,l_high_temp,500);
coeff= fgetl(f1);
coeff1= findstr(coeff,'E');
%extracting the coefficients
a1= coeff(1:coeff1(1)+3);
a1= str2num(a1);
a2= coeff(coeff1(1)+4:coeff1(2)+3);
a2= str2num(a2);
a3= coeff(coeff1(2)+4:coeff1(3)+3);
a3= str2num(a3);
a4= coeff(coeff1(3)+4:coeff1(4)+3);
a4= str2num(a4);
a5= coeff(coeff1(4)+4:coeff1(5)+3);
a5= str2num(a5);
new_line= fgetl(f1);
coeff2= findstr(new_line,'E');
a6= new_line(1:coeff2(1)+3);
a6= str2num(a6);
a7= new_line(coeff2(1)+4:coeff2(2)+3);
a7= str2num(a7);
a8= new_line(coeff2(2)+4:coeff2(3)+3);
a8= str2num(a8);
a9= new_line(coeff2(3)+4:coeff2(4)+3);
a9= str2num(a9);
a10= new_line(coeff2(4)+4:coeff2(5)+3);
a10= str2num(a10);
new_line1= fgetl(f1);
coeff_3= findstr(new_line1,'E');
a11= new_line1(1:coeff_3(1)+3);
a11= str2num(a11);
a12= new_line1(coeff_3(1)+4:coeff_3(2)+3);
a12= str2num(a12);
a13= new_line1(coeff_3(2)+4:coeff_3(3)+3);
a13= str2num(a13);
a14= new_line1(coeff_3(3)+4:coeff_3(4)+3);
a14= str2num(a14);
%specific heat calculation
Cp= specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,temp_range,R);
%enthalpy calculation
H= enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a2,a13,temp,temp_range,R);
%entropy calculation
S= entropy_1(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,temp_range,R);
%calculation of molecular weight
m_weight= weight(species_name);
%saving molecular weight
species_w= fopen('species_mass.txt','w');
fprintf(species_w,'molecular mass of %s is %f',species_name,m_weight);
fclose(species_w);
%plotting and saving
cur_dir=pwd;
%making a new folder for each species
mkdir(species_name);
cd(species_name);
figure(1)
plot(temp_range,Cp,'linewidth',2,'color',[0.2345 0.4567 0.5678])
xlabel('Temperature(K)');
ylabel('Specific heat(KJ/kg-K)');
title(species_name);
figure(2)
plot(temp_range,H,'LineWidth',2,'color',[0.5678 0.3423 0.7874])
xlabel('Temperature(K)');
ylabel('Enthalpy(KJ/Kg)');
title(species_name);
figure(3)
plot(temp_range,S,'LineWidth',2,'color',[0.4532 0.1234 0.3324])
xlabel('Temperature(K)');
ylabel('Entropy(KJ/K)');
title(species_name);
%saving the plots in spearate folder for each species
saveas(1,'specific heat.jpg');
saveas(2,'enthalpy.jpg');
saveas(3,'entropy.jpg');
cd(cur_dir);
%printing the molcular weight of the species
sprintf('molecular weight of %s = %f',species_name,m_weight)
end
%closing the file
fclose(f1)
WRITING THE FUNCTION FOR SPECIFIC HEAT
1. Now we will find the Specific heat by creating Specific heat function and assigning its resulted formula from governing equation namely, coefficent from 14 coefficients, Gas constant, Temperature and mid temperature.
2. Now we will create a condition using if loop and make the temperature less than and equal to mid temperature to specify lower temperature coefficent in the formula.
3. Then we will use else condition to specify coefficient of Higher temperature coefficent in the formula.
4. Now, we end the loop and also the function to calculate Specific Heat.
FUNCTION CODE FOR SPECIFIC HEAT:
%code for calculation of specific heat
function[Cp]= specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,temp_range,R)
%a1 a2 a3 a4 a5 are high temperature coefficients
%a8 a9 a10 a11 a12 are low temperature coefficients
for i= 1:length(temp_range)
if temp_range(i)<=temp(2)
Cp(i)= R*(a8 + (a9*temp_range(i)) + (a10*(temp_range(i).^2)) + (a11*(temp_range(i).^3)) + (a12*(temp_range(i).^4)));
else
Cp(i)= R*(a1 + (a2*temp_range(i)) + (a3*(temp_range(i).^2)) + (a4*(temp_range(i).^3)) + (a5*(temp_range(i).^4)));
end
end
WRITING THE FUNCTION FOR ENTHALPY
1. We will find the Enthalpy by creating Enthalpy function and assigning its relted formula coefficent from 14 coefficients, Gas constant, Temperature and mid temperature.
2. Now we will creatre an condition using if loop and make the temperature less than and equal to mid temperature to specify lower temperature coefficent in the formula.
3. Then we will use else condition to specify coefficient of Higher temperature coefficent in the formula.
4. Next, we end the loop and also the function to calculate Enthalpy.
FUNCTION CODE FOR ENTHALPY
%code for calculation of enthalpy
function[H]= enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,temp,temp_range,R)
%a1 a2 a3 a4 a5 are high temperature coefficients
%a8 a9 a10 a11 a12 are low temperature coefficients
for i= 1:length(temp_range)
if temp_range(i)<=temp(2)
H(i)= R*((a8*temp_range(i)) + (a9*((temp_range(i).^2)/2)) + (a10*((temp_range(i).^3)/3)) + (a11*((temp_range(i).^4)/4)) + (a12*((temp_range(i).^5)/5)) + a13);
else
H(i)= R*((a1*temp_range(i)) + (a2*((temp_range(i).^2)/2)) + (a3*((temp_range(i).^3)/3)) + (a4*((temp_range(i).^4)/4)) + (a5*((temp_range(i).^5)/5)) + a6);
end
end
WRITING THE FUNCTION FOR ENTROPY
1. Now we will find the Entropy by creating Entropy function and assigning its resulting formula coefficent from 14 coefficients, Gas constant, Temperature and mid temperature.
2. Now we will creatre an condition using if loop and make the temperature less than and equal to mid temperature to specify lower temperature coefficent in the formula.
3. Then we will use else condition to specify coefficient of Higher temperature coefficent in the formula.
4. Nextt, we will end the loop and also the function to calculate Entropy.
FUNCTION CODE FOR ENTROPY
%code for calculation of entropy
function[S]= entropy_1(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,temp_range,R)
%a1 a2 a3 a4 a5 are high temperature coefficients
%a8 a9 a10 a11 a12 are low temperature coefficients
for i= 1:length(temp_range)
if temp_range(i)<=temp(2)
S(i)= R.*((a8*log(temp_range(i))) + (a9*temp_range(i)) + (a10*((temp_range(i).^2)/2)) + (a11*((temp_range(i).^3)/3)) + (a12*((temp_range(i).^4)/4)) + a14);
else
S(i)= R.*((a1*log(temp_range(i))) + (a2*temp_range(i)) + (a3*((temp_range(i).^2)/2)) + (a4*((temp_range(i).^3)/3)) + (a5*((temp_range(i).^4)/4)) + a7);
end
end
WRITING THE FUNCTION FOR MOLECULAR WEIGHT
1. First, we create an Molecular_weight functionW and assigning species input to that.
2. Then define the main element of all 53 species in an array atomic_name that are :- 'H','C','O','N','Ar'.
3. Then by taking above main elements we will create there atomic_wieght array from any source like internet that are :- 1.008 12.011 15.999 14.007 39.948.Then we will assign molecular starting range W To zero and then create an for loop condition which runs to the length of the species and an another nested for loop is created which runs to the length of the string atomic_name.
4. Then we will create an condition using if condition and use strcmp command to compare the species and atomic_name strings to calculate molecular masses if by comparison both string found same then we will assign its molecular weight.
5. Now, we will string to the number using str2num command assign to n. Then we run an condition of if loop to find weather the species contains more than one element if there are more then 1 elements that is n>1 then we will increment the atomic weight by the multiplication of position and n-1 to the atomic_weight array.
6. Then we will print the Molecular weight using fprintf command and display it using disp command.
FUNCTION CODE FOR MOLECULAR WEIGHT
%code for calculation of molecular weight
function[m_weight]= weight(species_name)
element= {'H','O','N','Ar','C','S'};
atomic_mass= [1.008, 15.999, 14.007, 39.948, 12.011, 32.065];
m_weight=0;
for i= 1:length(species_name)
for j= 1:length(element)
if strcmp(species_name(i),element(j))
m_weight= m_weight + atomic_mass(j);
a= j;
end
end
b= str2num(species_name(i));
if b>1
m_weight= m_weight + atomic_mass(a)*(b-1);
end
end
endPLOTS SUPPORTING THE CODE
1. O2
2. N2
3. CO2
SCREENSHOT SHOWING MOLECULAR WEIGHT OF DIFFERENT SPECIES
Similarly, it goes on for 53 elements.
FOLDER showing 53 ELEMENTS containing Cp,S,E:
References:
1. https://en.wikipedia.org/wiki/Parsing
2. https://www.mathworks.com/matlabcentral/answers/317006-create-new-folders-by-mkdir
3. https://www.mathworks.com/matlabcentral/answers/106320-make-a-loop-to-calculate-molecular-weight
4. https://www.mathworks.com/matlabcentral/answers/17580-how-to-define-a-path-in-saveas-command
ERRORS OCCURED
Typing error resulted in unsuccessful run of the programme.
OVERCOMING THE ERROR: Mistake was rectified and the program did run successfully.