Project 1 - Parsing NASA thermodynamic data
Aim:
Write a function that extracts the 14 co-efficient and calculates the enthalpy, entropy and specific heats for all the species in the data file.
Objective:
- Calculate the molecular weight of each series
- Plot the Cp, enthalpy and entropy for the local temperature range (low temperature: high temperature) specific for each species.
Concept:
File Parsing:
It is the method of splitting the file or other input into a piece of data that can be easily stored or manipulated. By doing this we are going to split a string into parts then recognizing the parts to convert it into something simpler than a string.
Her in this project we are going to split the NASA thermodynamic data which is having 53 gas species value in 213 lines. To perform the calculation for thermodynamic properties, we need some sort of co-efficient values and temperature ranges which is available in the NASA file. So, we are going to parse the values which are all needed for our calculation from that file.
Specific heat:
Specific heat is the amount of heat energy required to raise the temperature of the body per unit mass. Specific heat is also known as specific heat capacity or mass specific heat. In SI units, specific heat is the amount of heat in joules required to raise 1 gm of a substance 1 kelvin.
Enthalpy:
When substance changes at constant pressure, enthalpy tells how much heat and work was added or removed from the substance. Enthalpy is similar to energy but not the same. When a substance grows or shrinks, energy is used up or released
Entropy:
The entropy of the object is a measure of the amount of energy which is unavailable to do work. Entropy is also measure of the number of possible arrangements of atoms a system can have. In this sense, entropy is the measure of uncertainty and randomness.
Gas Constant:
The gas constant is a physical constant denoted by R and is expressed in terms of units of energy per temperature increment per mole. The gas constant value is equivalent to Boltzmann constant but expressed as the pressure-volume produce instead of energy per increment of temperature per particle.
Molecular weight:
Molecular weight is a measure of sum of the atomic weight values of the atoms in a molecule. Molecular weight is used in chemistry to determine stoichiometry in chemical reactions and equations.
NASA thermodynamics file which is THERMO.dat file in which there are temperature range and 53 species are given in these data and each species contains 14 coefficients and temperature ranges and we have to parse these 14 coefficients to find the specific heat, entropy and enthalpy using formula as:
Where,
R = Universal gas constant
T = Temperature
Cp = Specific Heat (K)
S = Entropy (kJ/mole-K)
H = Enthalpy (kJ/mole)
a(i) = Coefficients where i 1 to 14
MATLAB Program:
clear all
close all
clc
%Inputs
%Gas Constant
R=8.314;
%Reading data file
%Open file and reading THERMO.dat file
f1=fopen('THERMO.dat','r');
get=fgetl(f1);
%Now reading temperature from this file
temp = fgetl(f1);
%splitting
s=strsplit(temp);
%converting strings cells into numbers
templow=str2double(s{2});
tempmid=str2double(s{3});
temphigh=str2double(s{4});
%Temperature range
T = linspace(templow,temphigh,1000);
%Skipping comment lines
for i=1:3
get=fgetl(f1);
end
%Reading the species
for j=1:53
tline=fgetl(f1);
A=strsplit(tline,' '); %splitting into cell
species = (A{1});
%To read the n
%Converting string to number
templow=str2double(A{end-3});
midtemp=str2double(A{end-2});
hightemp=str2double(A{end-1});
tlinel=fgetl(f1);
%getting next line
a=strfind(tlinel,'E');
%FInding first 7 higher and second 7 lower
%Finding Higher coefficients
a1=tlinel(1:a(1)+3);
a1=str2double(a1);
a2=tlinel(a(1)+4:a(2)+3);
a2=str2double(a2);
a3=tlinel(a(2)+4:a(3)+3);
a3=str2double(a3);
a4=tlinel(a(3)+4:a(4)+3);
a4=str2double(a4);
a5=tlinel(a(4)+4:a(5)+3);
a5=str2double(a5);
tline2=fgetl(f1);
b=strfind(tline2,'E');
a6=tline2(1:b(1)+3);
a6=str2double(a6);
a7=tline2(b(1)+4:b(2)+3);
a7=str2double(a7);
%Finding lower coefficients
a8=tline2(b(2)+4:b(3)+3);
a8=str2double(a8);
a9=tline2(b(3)+4:b(4)+3);
a9=str2double(a9);
a10=tline2(b(4)+4:b(5)+3);
a10=str2double(a10);
tline3 = fgetl(f1);
c=strfind(tline3,'E');
a11=tline3(1:c(1)+3);
a11=str2double(a11);
a12=tline3(c(1)+4:c(2)+3);
a12=str2double(a12);
a13=tline3(c(2)+4:c(3)+3);
a13=str2double(a13);
a14=tline3(c(3)+4:c(4)+3);
a14=str2double(a14);
%calculate specific heat
Cp=Specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,tempmid);
S=Entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,tempmid);
%calculate enthalpy
H=Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,R,tempmid);
%calculate molecular weight of the species
W=Molecular_weight(species);
%writing molecular weight file
MW_file = fopen('Molecular weight.txt','w');
fprintf(MW_file,'Molecular weight of %s is %d',species,W);
fclose(MW_file); %closing the file
%Plotting Cp, H and S with temperature
%Saving in folder
cd('C:\Users\udayn\Desktop\Matlab_challenges\Project_1'); %changing directory
mkdir(species); %creating directory
cd(species);
%Plot Cp vs T
figure(1)
plot(T,Cp,'linewidth',2,'color','b');
xlabel('Temperature(K)');
ylabel('Specific Heat(kJ)');
title(sprintf('Specific heat vs Temperature range of %s', species));
grid on
saveas (figure (1),'Specific heat.jpg');
%plot2 S vs T
figure(2)
plot(T,S,'linewidth',2,'color','r');
xlabel('Temperature');
ylabel('Entropy(kJ/mol-K)');
title(sprintf('Entropy vs Temperature range of %s',species));
grid on
saveas (figure(2),'Entropy.jpg');
%plot3 Cp vs T
figure (3)
plot(T,H,'linewidth',2,'color','g');
xlabel('Temperature');
ylabel('Enthalpy(kJ/mol)');
title(sprintf('Enthalpy vs Temperature range of %s',species))
grid on
saveas (figure(3),'Enthalpy.jpg');
cd('C:\Users\udayn\Desktop\Matlab_challenges\Project_1') %changing directory to original path
end
Code explanation:
- Use the input parameter R and we will open a file THERMO.dat file and reading this file using fopen command.
- Then we will read first THERMO using fgetl command and then we read temperature file and we will split the temperatures using strsplit command to convert it into cell array and then we will convert this string to number using str2double (str2num) command to find low, high, and medium temperatures.
- Then we will assign temperature range using linspace command. Then run loop to skip command lines the data file using fgetl command inside the loop.
- Now, we will run the loop to find 53 elements in the data file by assigning fgetl command and then splitting 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.
- Then we read the new line. In order to read the coefficient in the data file we will first find the E string using strfind command to find before and after numbers of the coefficients then we will find remaining 14 coefficients and convert them to numbers from strings. Out of which we will find first 7 coefficients which are for higher temperature and then second 7 coefficients which are for lower temperature.
- Now we will specify heat Cp, Entropy S, Enthalpy H function to find these parameters using these 14 coefficients for lower and higher temperature, then we will specify molecular weight W function to find the Molecular weight of the different species.
- Now we will write the Molecular weight.txt file and fprintf all the molecular weights of the species and then close the txt file using fclose command.
- Now we change the directory using cd command and then creating the directory folder where we have to save our files using mkdir command and then we will change the directory to this created folder. Now, we will plot specific heat, entropy and enthalpy with the temperature for different species
- Then we will use saveas command to save plots of different species in the directory folder. Now we willchnage the directory to which our main code file runs in the MATLAB.
Specific heat function code:
function Cp=Specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,tempmid)
%calculating specific heat
%claculating specific heat for higher and lower temoerature respectively
if (T <=tempmid)
Cp = R*(a1+(a2*T)+(a3*(T).^2)+(a4*(T).^3)+(a5*(T).^4));
else
Cp = R*(a8+(a9*T)+(a10*(T).^2)+(a11*(T).^3)+(a12*(T).^4));
end
- Now we will find the specific heat by creating specific heat function and assigning its related formula coefficient from 14 coefficients, Gas constant, Temperature and mid temperature.
- Now we will create a condition using if loop and make the temperature less than and equal to mind temperature to specify lower temperature coefficient in the formula. Then we will use else condition to specify coefficient of higher temperature coefficient in the formula.
- Then the end loop and function to calculate Specific Heat
Entropy
function S = Entropy (a1, a2, a3, a4, a5, a7, a8, a9, a10, a11, a12, a14, T, R, tempmid)
%calculating entropy
%calculating entropy for higher and lower temperature respectively
if (T <= tempmid)
S=R.*(a1.*log(T)+(a2.*T)+((a3*(T).^2/2)+((a4*(T).^3)/3)+((a5.*(T).^4)/4)+a7));
else
S = R.*(a8.*log(T)+(a9.*T)+((a10*(T).^2)/2)+((a11*(T).^3)/3)+((a12.*(T).^4)/4)+a14);
end
Now we will find the entropy by creating entropy function and assigning its related formula coefficient from 14 coefficients, Gas constant,
Temperature and mid temperature
- Now we will create a condition using if loop and make the temperature less than and equal to mid temperature to specify lower temperature coefficient in the formula.
- Then we will use else condition to specify coefficient of Higher temperature coefficient in the formula
- Then end the loop and function to calculate entropy.
Enthalpy:
function H = Enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,R,tempmid)
%calculating enthalpy
%calculating enthalpy for higher and lower temperature respectively
if(T<=tempmid)
H = R*T.*(a8+((a9*T)/2)+((a10*(T).^2)/3)+((a11*(T).^3)/4)+((a12*(T).^4)/5)+((a13./T)));
else
H = R*T.*(a1+((a2*T)/2)+((a3*(T).^2)/3)+((a4*(T).^3)/4)+((a5*(T).^4)/5)+((a6./T)));
end
end
Now we will find the Enthalpy by creating enthalpy function and assigning its related formula coefficient from 14 coefficients, Gas constant, temperature and mid temperature
- Now we will create a condition using if loop and make the temperature less than and equal to mid temperature to specify lower temperature coefficient in the formula.
- Then we will use else condition to specify coefficient of higher temperature coefficient in the formula.
- The end loop and function to calculate enthalpy
Molecular weight:
function W = Molecular_weight(species)
%Definig atomic name of the species
Atomic_name = ['H','C','O','N','Ar']
%Defining atomic masses of main species
atomic_weight = [1.008 12.011 15.999 14.007 39.948];
W=0; %Molecular weight starting range
for i=1:length(species) %loop runs to the length of the species
for j=1:length(Atomic_name) %loop runs to the length of the atomic names
if strcmp(species(i), Atomic_name(j)) %function to compare 2 strings
W = W+atomic_weight(j);
position = j;
end
end
n = str2double(species); %Now we will find if there is more element in the species then incrementing weight of the species
if n>1
W = W +atomic_weight(position)*(n-1);
end
end
fprintf('Molecular_weight of %s:-',species);
fprintf('%f',W);
disp('');
end
Firstly, we create a Molecular_weight function W and assigning species input to that.
- Then define the main element of all 53 species in an array atomic_name that are:- ‘HCO N Ar’
- Then by taking above main elements we will create there atomic_weight array from any source like internet that are:- 1.008 12.011 15.999 14.007 29.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.
- 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.
- Now we will string to the number using str2double 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 element than 1 that is n>1 then we will increment the atomic weight by the multiplication of position and n-1 to the atomic_weight array.
- Then we will print the Molecular weight using fprintf command and display it using disp command.
Command window output:
Atomic_name =
'HCONAr'
Molecular_weight of O:-15.999000
Atomic_name =
'HCONAr'
Molecular_weight of O2:-15.999000
Atomic_name =
'HCONAr'
Molecular_weight of H:-1.008000
Atomic_name =
'HCONAr'
Molecular_weight of H2:-1.008000
Atomic_name =
'HCONAr'
Molecular_weight of OH:-17.007000
Atomic_name =
'HCONAr'
Molecular_weight of H2O:-17.007000
Atomic_name =
'HCONAr'
Molecular_weight of HO2:-17.007000
Atomic_name =
'HCONAr'
Molecular_weight of H2O2:-17.007000
Atomic_name =
'HCONAr'
Molecular_weight of C:-12.011000
Atomic_name =
'HCONAr'
Molecular_weight of CH:-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of CH2:-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of CH2(S):-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of CH3:-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of CH4:-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of CO:-28.010000
Atomic_name =
'HCONAr'
Molecular_weight of CO2:-28.010000
Atomic_name =
'HCONAr'
Molecular_weight of HCO:-29.018000
Atomic_name =
'HCONAr'
Molecular_weight of CH2O:-29.018000
Atomic_name =
'HCONAr'
Molecular_weight of CH2OH:-30.026000
Atomic_name =
'HCONAr'
Molecular_weight of CH3O:-29.018000
Atomic_name =
'HCONAr'
Molecular_weight of CH3OH:-30.026000
Atomic_name =
'HCONAr'
Molecular_weight of C2H:-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of C2H2:-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of C2H3:-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of C2H4:-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of C2H5:-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of C2H6:-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of CH2CO:-41.029000
Atomic_name =
'HCONAr'
Molecular_weight of HCCO:-41.029000
Atomic_name =
'HCONAr'
Molecular_weight of HCCOH:-42.037000
Atomic_name =
'HCONAr'
Molecular_weight of H2CN:-27.026000
Atomic_name =
'HCONAr'
Molecular_weight of HCN:-27.026000
Atomic_name =
'HCONAr'
Molecular_weight of HNO:-31.014000
Atomic_name =
'HCONAr'
Molecular_weight of N:-14.007000
Atomic_name =
'HCONAr'
Molecular_weight of NNH:-29.022000
Atomic_name =
'HCONAr'
Molecular_weight of N2O:-30.006000
Atomic_name =
'HCONAr'
Molecular_weight of NH:-15.015000
Atomic_name =
'HCONAr'
Molecular_weight of NH2:-15.015000
Atomic_name =
'HCONAr'
Molecular_weight of NH3:-15.015000
Atomic_name =
'HCONAr'
Molecular_weight of NO:-30.006000
Atomic_name =
'HCONAr'
Molecular_weight of NO2:-30.006000
Atomic_name =
'HCONAr'
Molecular_weight of HCNO:-43.025000
Atomic_name =
'HCONAr'
Molecular_weight of HOCN:-43.025000
Atomic_name =
'HCONAr'
Molecular_weight of HNCO:-43.025000
Atomic_name =
'HCONAr'
Molecular_weight of NCO:-42.017000
Atomic_name =
'HCONAr'
Molecular_weight of CN:-26.018000
Atomic_name =
'HCONAr'
Molecular_weight of HCNN:-41.033000
Atomic_name =
'HCONAr'
Molecular_weight of N2:-14.007000
Atomic_name =
'HCONAr'
Molecular_weight of AR:-39.948000
Atomic_name =
'HCONAr'
Molecular_weight of C3H8:-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of C3H7:-13.019000
Atomic_name =
'HCONAr'
Molecular_weight of CH3CHO:-42.037000
Atomic_name =
'HCONAr'
Molecular_weight of CH2CHO:-42.037000>>
Plots of some species:
Species O2 plot enthalpy vs temp
Plots O2 entropy vs temp:
Plots species O2 specific heat vs temp
Species CO2 plot Enthalpy vs temp:
Species CO2 plot Entropy vs Temp
Species CO2 plot specific heat vs temp
Species N2 plot Enthalpy vs Temp
Species N2 plot Entropy vs Temp
Species N2 plot Specific Heat vs Temp
Graph file of all species
