Project 1 - Parsing NASA thermodynamic data
Aim:-
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. You will be reading this file NASA thermodynamic data(Tip: use fopen and fgetl). The formulae are shown below.
Here R is the universal gas constant, T is the temeprature
The FIRST 7 coefficients are HIGH-temperature coefficients and the SECOND 7 coefficients are LOW-temperature coefficients.
Also use the LOCAL TEMPERATURE RANGES which is unique for each individual species, not the GLOBAL TEMPERATURE RANGE.
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.
6. Please write a detailed report to describe your logic.
Main code:-
close all
clear all
clc
f1 = fopen('THERMO.dat','r');
fgetl(f1);
% Defining temperature values
templ = str2num(fgetl(f1)); %#ok<ST2NM>
global_low_temp = templ(1);
global_mid_temp = templ(2);
global_high_temp = templ(3);
% For skipping unwanted lines
for i = 1:3
fgetl(f1);
end
% Loop creation to read the values and plot molecular weight
for i = 1:4:212
line1 = fgetl(f1);
l = strsplit(line1, ' ');
m = length(l);
species = l{1};
% Storing the values of the local temperature values
low_temp = str2double(l(m-3));
mid_temp = str2double(l(m-2));
high_temp = str2double(l(m-1)); % Check this later
% Storing the values of the second line of variables
line2 = fgetl(f1);
% Finding the position of E
a = strfind(line2,'E');
% Reading all the other variables
a1 = str2double(line2(1:(a(1)+3)));
a2 = str2double(line2((a(1)+4):(a(2)+3)));
a3 = str2double(line2((a(2)+4):(a(3)+3)));
a4 = str2double(line2((a(3)+4):(a(4)+3)));
a5 = str2double(line2((a(3)+4):(a(5)+3)));
% Storing the values of the 3rd line
line3 = fgetl(f1);
% Finding the position of E
b = strfind(line3,'E');
% Storing the values for the next set of variables
a6 = str2double(line3(1:b(1)+3));
a7 = str2double(line3((b(1)+4):(b(2)+3)));
a8 = str2double(line3((b(2)+4):(b(3)+3)));
a9 = str2double(line3((b(3)+4):(b(4)+3)));
a10 = str2double(line3((b(4)+4):(b(5)+3)));
% Storing ther values of the 4th line
line4 = fgetl(f1);
% Finding the position of E
c = strfind(line4,'E');
% Storing the values of variables in line 4
a11 = str2double(line4(1:c(1)+3));
a12 = str2double(line4((c(1)+4):(c(2)+3)));
a13 = str2double(line4((c(2)+4):(c(3)+3)));
a14 = str2double(line4((c(2)+4):(c(4)+3)));
% Creating Temperature values for lower and higher temperature
lower_temp = linspace(low_temp,mid_temp,100);
higher_temp = linspace(mid_temp,high_temp,100);
% Value of the universal has constant in j/mol*K
R = 8.314; % Universal Gas Constant
% Calculating the values for cp, H and S using the given formulas
for k = 1:100
% lower temperature values for cp,h,s calculations
cp_low(k)=(a8+(a9*lower_temp(k))+(a10*(lower_temp(k)^2)+(a11*(lower_temp(k))^3)+a12*(lower_temp(k))^4))*R;
h_low(k)=(a8+(a9*lower_temp(k))/2+(a10*lower_temp(k)^2)/3+(a11*lower_temp(k)^3)/4+a12*(lower_temp(k)^4)/5+(a13/lower_temp(k)))*R*lower_temp(k);
s_low(k)=((a8*log(lower_temp(k))+(a9*lower_temp(k))+((a10*lower_temp(k)^2)/2)+((a11*lower_temp(k)^3)/3)+((a12*lower_temp(k)^4)/4)+a14)*R);
% higher temperature values for cp,h,s calculations
cp_high(k)=(a1+(a2*higher_temp(k))+(a3*(higher_temp(k)^2)+(a4*(higher_temp(k))^3)+a5*(higher_temp(k))^4))*R;
h_high(k)=(a1+(a2*higher_temp(k))/2+(a3*higher_temp(k)^2)/3+(a4*higher_temp(k)^3)/4+(a5*(higher_temp(k)^4))/5+(a6/higher_temp(k)))*R*higher_temp(k);
s_high(k)=((a1*log(higher_temp(k))+a2*higher_temp(k)+((a3*higher_temp(k)^2)/2)+((a4*higher_temp(k)^3)/3)+((a5*higher_temp(k)^4)/4)+a7)*R);
end
% Making Folders for saving the plots
nasa_graphs = species;
mkdir(nasa_graphs);
current = pwd;
cd(nasa_graphs);
% Plotting the graphs
% For specific heat
figure(1)
plot(lower_temp,cp_low,'linewidth',3)
hold on
plot(higher_temp,cp_high,'linewidth',3)
hold off
title(['Temperature vs Specific Heat',species])
xlabel('Temperature');
ylabel('Specific Heat');
legend('low temperature','high temperature');
fig = sprintf('Specific Heat of %s.png',species);
saveas(gca,fig);
% For Enthalpy
figure(2)
plot(lower_temp,h_low,'linewidth',3)
hold on
plot(higher_temp,h_high,'linewidth',3)
hold off
title(['Temperature vs Enthalpy',species])
xlabel('Temperature');
ylabel('Enthalpy');
legend('low temperature','high temperature');
fig = sprintf('Enthalpy of %s.png',species);
saveas(gca,fig);
% For Entropy
figure(3)
plot(lower_temp,s_low,'linewidth',2)
hold on
plot(higher_temp,s_high,'linewidth',2)
hold off
title(['Temperature vs Entropy',species])
xlabel('Temperature');
ylabel('Entropy');
legend('low temperature','high temperature');
fig = sprintf('Entropy of %s.png',species);
saveas(gca,fig);
cd(current)
end
% Calculating molecular weight
atoms = ['H','C','N','R','S'];
atom_weight = [1.00784 12.0109 16.002 39.065 32.067];
mol_weight = 0;
species_check = char(species);
for i=1:length(species_check)
for j=1:length(atoms)
if (strcmp(species_check(i),atoms(j)))
mol_weight = mol_weight + atom_weight(j);
position = j;
end
end
x = str2double(species_check(i));
if x>1
mol_weight = mol_weight + (atom_weight(j)*(x-1));
end
fprintf('Molecular weight of %s is : %f \n',species_check,mol_weight)
end EXPLAINATION:
1. The code uses the thermodynamic data values from the file thermo.dat that was provided to me.
2. This file was read with the help of the read command for files and was read line by line.
3. For getting values for various species the program used loops for manipulating the values to an extent of the number of iterations that have to be taken into account.
4. the co-efficients are arranged ina fixed order and we have found the exact positiions with the help of strfind command in matlab and using the find variable to be E.
5. We have made a folder with the help of mkdir command
6. We have assigned certain values to the variables in the file and plotted themp against various temperature values vs enthalpy, entropy.
7. The images were saved in .JPEG format with the help of the saveas command.
PLOTS:
The first set of plots include the plots for CO2.
TEMPERATURE VS ENTHALPY
TEMPERATURE VS SPECIFIC HEAT
TEMPERATURE VS ENTROPY