%Program No - 1 -writing a Program to perform the parsing of NASA Thermodynamic data of 53
%species ,to extract the co-efficients of the species and to find the
%thermodynamic properties such as Specific heat Capacity(Cp),Enthalpy(H)
%and Entropy(S).To plot the graphs and also to find the molecular weight of
%each specie.
close all
clear all
clc
%Opening the file which contains the Thermodynamic data
TD_data = fopen('THERMO data.txt','r') ; %open the file which is in txt format in read mode
name_file_header = fgetl(TD_data); % getting the first line of the file which is the header
Temperature_line= fgetl(TD_data); %getting the next line which contains the global temperature range
Global_Temperature_range = strsplit(Temperature_line,' '); %splitting the line to extract the global temperature numbers in the form of a string
Global_low_Temperature = str2num(Global_Temperature_range{2}) ; %converting the strings into numbers
Global_mid_Temperature = str2num(Global_Temperature_range{3}) ;
Global_high_Temperature = str2num(Global_Temperature_range{4}) ;
%assigning the numbers to the global temperature variables
GTl = Global_low_Temperature ; % low
GTm = Global_mid_Temperature ;%mid
GTh = Global_high_Temperature ;%high
%extracting the sentences
skip_lines1 = fgetl(TD_data);
skip_lines2 = fgetl(TD_data);
skip_lines3 = fgetl(TD_data);
%for-loop is used to extract the co-efficients of each of 53 species and
%find the Cp ,H ,S values as well plot them against local temperature
%range.And also find the molecular weight
for i = 1:53
introline = fgetl(TD_data); %getting the species represented line which has the details such as local temperature range
Local_Temp = strsplit(introline,' '); %splitting the line to extract the local temperatures
%using for loop to extract the local temperature values - Tl = low temperature,Th = High temperature,Tm = mid temperature every
%time and converting them into numbers , locating them in
%reverse direction and converting to numbers
for n=length(Local_Temp)
elemental_specie = Local_Temp{1} ; % the first string character is the name of specie
Tl = str2num(Local_Temp{n-3}); %low
Th = str2num(Local_Temp{n-2}); %high
Tm = str2num(Local_Temp{n-1}); %mid
end
R= 8.3145 ; % Universal gas constant units - J/mol-K
T = linspace(Tl,Th,500); %the local temperature in the form of an array , to find the property values as well as plot the graphs
%Extracting the polynomial co-efficients which are spread across next three
%lines in the data file where first two lines have 5 co-efficients and the
%last one has 4 coefficients. Getting these lines and then extracting them
%to convert into numbers
cf_line_1=fgetl(TD_data);
cf1= strfind(cf_line_1,'E'); %locating E to differentiate between the between the co-efficients
k1= str2num(cf_line_1(1:cf1(1)+3));%finding and converting them to numbers
k2= str2num(cf_line_1(cf1(1)+4:cf1(2)+3));
k3= str2num(cf_line_1(cf1(2)+4:cf1(3)+3));
k4= str2num(cf_line_1(cf1(3)+4:cf1(4)+3));
k5= str2num(cf_line_1(cf1(4)+4:cf1(5)+3));
cf_line_2=fgetl(TD_data);
cf2= strfind(cf_line_2,'E');
k6 = str2num(cf_line_2(1:cf2(1)+3));
k7 = str2num(cf_line_2(cf2(1)+4:cf2(2)+3));
k8 = str2num(cf_line_2(cf2(2)+4:cf2(3)+3));
k9 = str2num(cf_line_2(cf2(3)+4:cf2(4)+3));
k10= str2num(cf_line_2(cf2(4)+4:cf2(5)+3));
cf_line_3=fgetl(TD_data);
cf3= strfind(cf_line_3,'E');
k11= str2num(cf_line_3(1:cf3(1)+3));
k12= str2num(cf_line_3(cf3(1)+4:cf3(2)+3));
k13= str2num(cf_line_3(cf3(2)+4:cf3(3)+3));
k14= str2num(cf_line_3(cf3(3)+4:cf3(4)+3));
%calling the function which calculates the values of thermodynamic
%properties such as Cp ,H and S which are of the polynomial functions and
%depends upon the co-efficients extracted above and the local temperature
%values. These are substituted in the expressions and the respective values
%are found
[Cp,H,S] = TD_ppt(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,T,Tm,R);
%creating a folder where the plots between temperature and thermodynamic
%properties are to be saved
current_folder = pwd ;
mkdir(elemental_specie);
cd(elemental_specie);
%Plotting
%plot between Temperature Vs Specific Heat Capacity
figure(1)
plot(T,Cp,'linewidth',3,'color','r')
xlabel('Temperature[K]')
ylabel('Specific Heat Capacity[J/mol-K]')
title(elemental_specie)
grid on
saveas(figure(1),'Temperature Vs Specific Heat Capacity.png');
%plot between Temperature Vs Enthalpy
figure(2)
plot(T,H,'linewidth',3,'color','b')
xlabel('Temperature[K]')
ylabel('Enthalpy[J/kg]')
title(elemental_specie)
grid on
saveas(figure(2),'Temperature Vs Enthalpy.png');
%plot between Temperature Vs Entropy
figure(3)
plot(T,S,'linewidth',3,'color','y')
xlabel('Temperature[K]')
ylabel('Entropy[J/K]')
title(elemental_specie)
grid on
saveas(figure(3),'Temperature Vs Entropy.png');
cd(current_folder);
%To calculate the molecular weights of the species,These elemental species
%are molecules which are made up of Hydrogen (H) , Carbon (C) ,Oxygen(O)
%,Nitrogen(N) ,Argon(AR).
elements = ['H','C','O','N','AR'] ;
atomic_weights = [1 12 16 14 40] ; %the atomic weights of the elements involved are mentioned respectively
molecular_weight=0 ;
%creating a for-loop to find the molecular weights and comparing the
%strings of elemental species, locating them and to check whether they match and to
%consider them . Use it in the expression to calculate the molecular
%weight by using its corresponding atomic weight
for j= 1:length(elemental_specie)
for k = 1:length(atomic_weights)
if strcmp(elemental_specie(j),elements(k))
molecular_weight= molecular_weight+atomic_weights(k) ;
location = k;
end
end
l = str2double(elemental_specie(j));
if l>1
molecular_weight=molecular_weight+atomic_weights(location)*(l-1) ;
end
end
fprintf('molecular_weight of %s - ',elemental_specie) %to print the sentence and display the value of the calculated molecular weight
fprintf('%d atomic mass units ',molecular_weight) %to display the units of molecular weight after the number
end
fclose(TD_data); %close the file and the program ends after this
%Program No - 2 -Writing a separate program for defining the function to find the
%thermodynamic properties of Cp,H and S .There are two sets of coefficients
%having 7 in each. The first 7 are to be used for higher temperatures in the
%properties expression while the next 7 are to be used for expressions
%dealing with lower temperatures.Tm is the local mid temperature and T= local
%temperature range from lower to higher. The properties are polynomials which depend upon local temperature
%values, co-efficients extracted from the Thermo polynomial data of NASA
% and Universal gas constant
% (R).(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14) = polynomial
% co-efficients
function [Cp,H,S] = TD_ppt(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,T,Tm,R)
for n = 1:length(T)
if T(n) >= Tm %using if clause statement to bifurcate using the mid temperature as extracted from the thermo data file
%the set of polynomial expressions corresponding to local high temperature
%range
Cp = R*(k1 + k2.*(T)+k3.*((T).^2)+k4.*((T).^3)+k5.*((T).^4)) ;
H = R.*(T).*(k1+(k2.*(T./2))+(k3.*(T.^2./3))+(k4.*(T.^3./4))+(k5.*(T.^4./5))+(k6./T)) ;
S = R*((k1.*log(T))+k2.*(T)+(k3.*(T.^2./2))+(k4.*(T.^3./3))+(k5.*(T.^4./4))+k7);
else
%the set of polynomial expressions corresponding to local low temperature
%range
Cp = R.*(k8 + k9.*(T)+k10.*((T).^2)+k11.*((T).^3)+k12.*((T).^4));
H = R.*(T).*(k8+(k9.*(T./2))+(k10.*(T.^2./3))+(k11.*(T.^3./4))+(k12.*(T.^4./5))+(k13./T)) ;
S = R*((k8.*log(T))+ k9.*(T)+(k10.*(T.^2./2))+(k11.*(T.^3./3))+(k12.*(T.^4./4))+k14) ;
end
end
end