To write a project to evalute the NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S. They have also documented the co-efficients that are required to evaluate these polynomials for 1000+ species.
1. 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. Calculate the molecular weight of each species and display it in the command window.
3. Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species.
4. 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.
5. The main program that you will submit will just take the "THERMO.dat" file as an input and generate Cp, H and S plots.
clc
% For getting the data file
f1=fopen('THERMO.dat','r');
fgetl(f1);
% Reading the first line of the file
header_global_temp=fgetl(f1);
% Splitting the string
A=strsplit(header_global_temp,' ');
% Obtaining global max,medium and minimum temperature:
global_min=strsplit(A{2});
global_mid=strsplit(A{3});
global_max=strsplit(A{4});
% ignoring the comments in the file
for i=1:3
fgetl(f1);
end
% Reading all line data
for i=1:53
header=fgetl(f1);
species_header=strsplit(header,' ');
species_name=species_header{1};
% Finding the string
loc_temp=strsplit(header,' ');
speciess_name=(loc_temp(1));
loc_min_temp=str2double(loc_temp(end-3));
loc_med_temp=str2double(loc_temp(end-2));
loc_max_temp=str2double(loc_temp(end-1));
% Finding the coefficients a1 to a14
tline1=fgetl(f1);
tline2=fgetl(f1);
tline3=fgetl(f1);
a=strfind(tline1,'E');
a1=tline1(1:a(1)+3);
a1=str2double(a1);
a2=tline1(a(1)+4:a(2)+3);
a2=str2double(a2);
a3=tline1(a(2)+4:a(3)+3);
a3=str2double(a3);
a4=tline1(a(3)+4:a(4)+3);
a4=str2double(a4);
a5=tline1(a(4)+4:a(5)+3);
a5=str2double(a5);
% Finding strings and convert to numbers
a=strfind(tline2,'E');
a6=tline2(1:a(1)+3);
a6=str2double(a6);
a7=tline2(a(1)+4:a(2)+3);
a7=str2double(a7);
a8=tline2(a(2)+4:a(3)+3);
a8=str2double(a8);
a9=tline2(a(3)+4:a(4)+3);
a9=str2double(a9);
a10=tline2(a(4)+4:a(5)+3);
a10=str2double(a10);
% Finding strings and convert to numbers
a=strfind(tline3,'E');
a11=tline3(1:a(1)+3);
a11=str2double(a11);
a12=tline3(a(1)+4:a(2)+3);
a12=str2double(a12);
a13=tline3(a(2)+4:a(3)+3);
a13=str2double(a13);
a14=tline3(a(3)+4:a(4)+3);
a14=str2double(a14);
% Creating folder
cd(['C:\Users\kavinish\OneDrive\Documents\MATLAB\',species_name]);
T=linspace(loc_min_temp,loc_max_temp,1000);
R=8.314;
% For getting the array of cp,h.and s
for k=1:length(T)
% Temparature to assign the coefficients
% a1-a7 is high temperature coefficients
% a8-a14 is low temperature coefficients
if T(k)>loc_med_temp
% high temparature coefficients
cp=(a1+a2.*T+a3.*T.^2+a4.*T.^3+a5.*T.^4)*R;
H=(a1+(a2.*T)/2+(a3.*T.^2)/3+(a4.*T.^3)/4+(a5.*T.^4)/5+(a6./T))*R.*T;
S=a1.*log(T)+a2.*T+(a3.*T.^2)/2+(a4.*T.^3)/3+(a5.*T.^4)/4+a7;
else
%low temparature coefficients
cp=(a8+a9.*T+a10.*T.^2+a11.*T.^3+a12.*T.^4)*R;
H=(a8+(a9.*T)/2+(a10.*T.^2)/3+(a11.*T.^3)/4+(a12.*T.^4)/5+(a13./T))*R.*T;
S=a8.*log(T)+a9.*T+(a10.*T.^2)/2+(a11.*T.^3)/3+(a12.*T.^4)/4+a14;
end
end
%for i=1:53
%loc_temp=strsplit(tline1,' ');
%species_name=(loc_temp(1));
% plotting the curves
%cp v/s T
figure(1);
plot(cp,T,'LineWidth',3,'Color','g');
grid on
title(species_name);
xlabel('specific heat');
ylabel('Temparature');
saveas(1,'cp_vs_temp.jpeg');
% H v/s T
figure(2);
plot(H,T,'LineWidth',3,'Color','r');
title(species_name);
xlabel('Enthalpy');
ylabel('Temparature');
% entropy v/s temparature
figure(3);
plot(S,T,'LineWidth',3,'Color','b');
title(species_name);
xlabel('Entropy');
ylabel('Temparature');
% changing the directory and saving the file
saving=(['C:\Users\kavinish\OneDrive\Documents\MATLAB\',species_name]);
cd(saving)
saveas(1,'cp_vs_temp.jpeg');
saveas(2,'enthalpy_vs_temp.jpeg');
saveas(3,'entropy_vs_temp.jpeg');
end
COMMAND WINDOW:
function[w]=mol_weight(species_name)
atomic_names=['H','C','O','N','A'];
atomic_weights=[1,12,16,14,40];
w=0;
for i=1:length(species_name)
for j=1:length(atomic_names)
if strcmpi(species_name(i),atomic_names(j))
w=w+atomic_weights(j);
pos=j;
end
n=str2num(species_name(i));
end
if n>=1
w=w+atomic_weights(pos)*(n-1);
k=fprintf('Molecular weight of %s = %f',species_name,w);
end
end