Parsing NASA thermodynamic data
AIM:
- 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.
- To Calculate the molecular weight of each species and to display it in the command window.
- To Plot the Cp, Enthalpy and Entropy for the local temperature range specific for each species.
- To Save the plots as images with appropriate names and dump them in separate folders for each species with suitable name.
DESCRIPTION:-
1.) Initially open the given thermodynamic data file and read the data in the file and store the data to the variable F.
2.) Then read the header file of thermodynamic data and calculate the global low,mid and high temperature from the second line of header file.
3.) By using fgetl(), strsplit() and str2num() command find out the global low,mid and high temperatures in the second line.
4.) Then calculate the 14 coefficients for all the species in the given thermodynamic data file.
5.) In data file each species has three lines in which 14 coefficients are there.
6.) By using the for loop, calculate the 14 coefficients for all the species.
7.) By using strfind() command, find out the number of exponential in first line, by which we can find out the number of coefficient in first line.
8.) Once you find out, then convert the coefficient from string form to number form by using str2num() command.
9.) By doing the above steps for line two and three, we can calculate the 14 coefficients for all the species.
10.) Once you find out the coefficients for a species, then calculate the entropy,enthalpy and specific heat.
11.) First 7 coefficients are high temperature coefficients and calculate cp,H and s.
12.) Assuming universal gas constant R=8.314 J/mol-k and assuming temperature range from 100 to 5000 take 100 values by linspace() command.
13.) Once you calculate the entropy,enthalpy and specific heat for all the species, plot the graph between: Temperature vs Entropy, Temperature vs Enthalpy and Temperature vs Specific heat.
14.) Store all the three plots for all the species in the file. The stored file is shown in the above picture.
15.) Finally calculating the molecular weight for all the species and displayed in the command window.
16.) The method for calculating molecular weight is shown in the function program.
GOVERNING EQUATIONS:-
CpR=a1+a2T+a3T2+a4T3+a5T4CpR=a1+a2T+a3T2+a4T3+a5T4
HRT=a1+a2T2+a3T23+a4T34+a5T45+a6THRT=a1+a2T2+a3T23+a4T34+a5T45+a6T
SR=a1lnT+a2T+a3T22+a4T33+a5T44+a7SR=a1lnT+a2T+a3T22+a4T33+a5T44+a7
CODES:-
funtion in code
function [M] = molecular_mass(species_name)
species = upper(species_name);
atoms = 'HCNO';
atomic_weight = [1.00784 12.0096 14.00643 15.99903];
M = 0;
for i=1:length(species)
n=str2num(species(i));
if n>1
M=M+p*(n-1);
else
for j=1:4
if strcmp(species(i), atoms(j))
M = M +atomic_weight(j);
p=atomic_weight(j);
end
end
end
end
for l = 1:length(species)
if strcmp(species(l),'R')
M = M + 38.8;
end
end
end
MAIN CODE:-
close all
clear all
clc
%open the NASA thermodynamic data in reading form
F = fopen('THERMO.dat','r')
%header file information
l1=fgetl(F);
l2=fgetl(F);
l3=fgetl(F);
l4=fgetl(F);
l5=fgetl(F);
%calculating low,mid,high global temperaure in line2(i.e.l2)
a=strsplit(l2,' ');
global_low_temp = str2num(a{2})
global_mid_temp = str2num(a{3})
global_high_temp = str2num(a{4})
%extracting 14 coefficient for each species
for i=1:53 %total no.of species = 53
Species= fgetl(F);
species_name = Species(1:6)
first_line = fgetl(F);
E1 = strfind(first_line,'E');
c1 = str2num(first_line(1:(E1(1)+3)));
c2 = str2num(first_line((E1(1)+4):(E1(2)+3)));
c3 = str2num(first_line((E1(2)+4):(E1(3)+3)));
c4 = str2num(first_line((E1(3)+4):(E1(4)+3)));
c5 = str2num(first_line((E1(4)+4):(E1(5)+3)));
second_line = fgetl(F);
E2 = strfind(second_line,'E');
c6 = str2num(second_line(1:(E2(1)+3)));
c7 = str2num(second_line((E2(1)+4):(E2(2)+3)));
c8 = str2num(second_line((E2(2)+4):(E2(3)+3)));
c9 = str2num(second_line((E2(3)+4):(E2(4)+3)));
c10 = str2num(second_line((E2(4)+4):(E2(5)+3)));
third_line = fgetl(F);
E3 = strfind(third_line,'E');
c11 = str2num(third_line(1:(E3(1)+3)));
c12 = str2num(third_line((E3(1)+4):(E3(2)+3)));
c13 = str2num(third_line((E3(2)+4):(E3(3)+3)));
c14 = str2num(third_line((E3(3)+4):(E3(4)+3)));
%creating folder for each species
mkdir(['C:Userssubramania sivaDesktopmatlab programsthermospecies',species_name]);
%assuming temperature ranges
temp = linspace(100,5000,100);
%universal gas constant
R=8.314;
%calculation of Cp,H and S for each species
for k=1:length(temp)
if temp(k)<3000
Cp = R*(c1+c2.*temp+c3.*temp.^2+c4.*temp.^3+c5.*temp.^4);
H = R.*temp.*(c1+(c2.*temp)./2+(c3.*temp.^2)/3+(c4.*temp.^3)/4+(c5.*temp.^4/5)+c6./temp);
S = R.*(c1*log(temp)+c2.*temp+(c3.*temp.^2)/2+(c4.*temp.^3)/3+(c5.*temp.^4)/4+c7);
else
Cp = R.*(c8+c9.*temp+c10.*temp.^2+c11.*temp.^3+c12.*temp.^4);
H = R.*temp.*(c8+(c9.*temp)/2+(c10.*temp.^2)/3+(c11.*temp.^3)/4+(c12.*temp.^4/5)+c13./temp);
S = R.*(c8.*log(temp)+c9.*temp+(c10.*temp.^2)/2+(c11.*temp.^3)/3+(c12.*temp.^4)/4+c14);
end
end
%plotting
%temperature vs entropy
figure(1);
plot(temp,S,'linewidth',3,'color','r');
title(species_name);
xlabel('Temperature');
ylabel('Entropy');
%temperature vs specific heat
figure(2);
plot(temp,Cp,'linewidth',3,'color','b');
title(species_name);
xlabel('Temperature');
ylabel('Specific Heat');
%Temperature vs Enthalpy
figure(3);
plot(temp,H,'linewidth',3,'color','g');
title(species_name);
xlabel('Temperature');
ylabel('Enthalpy');
%changing the directory to save the plots
address =(['C:Userssubramania sivaDesktopmatlab programsthermospecies',species_name]);
cd(address);
%saving the plots to its respective species
saveas(1,'Entropy.jpg');
saveas(2,'specific heat.jpg');
saveas(3,'Enthalpy.jpg');
%changing the directory for calculating the molecular_weight
cd('C:Userssubramania sivaDesktopmatlab programsthermo')
Molecular_weight = molecular_mass(species_name)
end
fclose(F)
OUTPUT:-
References: