Project 1 - Parsing NASA thermodynamic data
Specific Heat: The heat required to raise the temperature of 1 gram of a substance to 1 degree Celsius.
we have specific heat at constant pressure(Cp) and specific heat at constant volume(Cv)
In order to code the Cp function we will be using the above-mentioned formula
%%making the function for specific heat calculation
%%function definition with the co-effiecients
function CP =s_h(c1,c2,c3,c4,c5,c8,c9,c10,c11,c12,T,R,temp_mid)
if (T>temp_mid)
CP = R*(c1+ (c2*T) + (c3*(T).^2) + (c4*(T).^3) + (c5*(T).^4));
else
CP =R*(c8+ (c9*T) + (c10*(T).^2) + (c11*(T).^3) + (c12*(T).^4));
end
endIn the above code, we will check the temperature range, and if the temperatures are high, we will use the first seven coefficients. We will consider the coefficients from 8-14 if we have low temperatures.
Enthalpy:
Enthalpy is the partial Legendre transform of U(internal energy) that replaces the volume by the pressure as the independent variable.
H = U(p)
H = U +pv
%Enthalpy calculation function
function H = Enthalpy(c1,c2,c3,c4,c5,c6,c8,c9,c10,c11,c12,c13,T,R,temp_mid) % function definition
if (T>temp_mid) % if condition to check the temperatutre
H = R*(T).*(c1 + ((c2*(T))/2) + ((c3*((T)).^2)/3) + ((c4*((T)).^3)/4) + ((c5*((T)).^4)/5) + ((c6./(T)))); % equation to calculate the enthalpy
else % if IF condition not satisfied
H =R*(T).*(c8 + ((c9*(T))/2) + ((c10*((T)).^2)/3) + ((c11*((T)).^3)/4) + ((c12*((T)).^4)/5) + ((c13./(T)))); % equation to calculate for low twmperatures <= temp_mid
end % end if
end % end functionIn the above code, we will check the temperature range, and if the temperatures are high, we will use the first seven coefficients. We will consider the coefficients from 8-14 if we have low temperatures.
Entropy: Entropy is a function of the extensive parameters of any composite system, defined for all equilibrium states and having the values assumed by extensive parameters in the absence of internal constraint.
% Entropy Calculation Function
function S = Entropy(c1,c2,c3,c4,c5,c6,c8,c9,c10,c11,c12,c13,T,R,temp_mid) % function definition
% calculating entropy for higher and lower temperature respectively.
if (T >temp_mid) % if condition to check the temperatutre
S = R*(c1.*log(T) + (c2*T) + ((c3*(T).^2)/2) + ((c4*(T).^3)/3) + ((c5*(T).^4)/4) +((c6*(T).^5)/5)) ; % equation to calculate the entropy
else % if IF condition not satisfied
S = R*(c8.*log(T) + (c9*T) + ((c10*(T).^2)/2) + ((c11*(T).^3)/3) + ((c12*(T).^4)/4) +((c13*(T).^5)/5)); % equation to calculate for low twmperatures <= temp_mid
end % end if
end % end functionIn the above code, we will check the temperature range, and if the temperatures are high, we will use the first seven coefficients. We will consider the coefficients from 8-14 if we have low temperatures.
Molecular weight_calculation:
In order to calculate the molecular weight we will use the following formula
MW = molar mass * number of moles of the element
% Molecular Weight Calculation Function
function MW = Molecular_weight(elements)
% Defining element Name
% Hydrogen , Carbon, Oxygen, Nitrogen, Argon
element_name = ['H', 'C', 'O', 'N', 'Ar'];
% Defining Atomic mass of element
atomic_mass = [1.008 12.011 15.999 14.007 39.948];
% Molecular weight initialvalue
MW = 0;
% Adding Atomic Number to obtain particular value
for i=1:length(elements) % Loop runs to the length of the elements
for j=1:length(element_name) % Loop runs to the length of the atomic name
if strcmp(elements(i), element_name(j)) % string comparision function to compare 2 strings
MW = MW + atomic_mass(j);
state = j;
end % end if
end % end for
% Finding more elements in the substance and increasing atomic weight of the elements
n = str2double(elements);
if n>1 % if condition if the number of moles are greater than 1
MW = MW + atomic_mass(state)*(n-1);
end % end if
end % end for
%printing the molecular weight of the element
fprintf('Molecular Weight of %s : ',elements);
fprintf('%f',MW);
disp(' ');
end % end functionThis function is used to calculate the molecular weight of each substance.
Main Code:
clear
close all
clc
% universal Gas Constant (J/(mol-K))
R = 8.314;
% Open and read THERMO.dat file
% Skip the first line (header)
f1=fopen('THERMO.dat','r');
fgetl(f1);
% reading the header information
termline =fgetl(f1);
% Splitting the temperature range line into 3 temperatures and asssigning
% the value
A =str2num(termline);
temp_low = A(:,1);
temp_mid = A(:,2);
temp_high = A(:,3);
%Temperature range for the values where we caliculate Cp,H,S
T =linspace(temp_low,temp_high,1000);
%Eleminating the commented lines
for i =1:3
fgetl(f1);
end
%Read the data of elements
for i = 1:53 % all 53 elements or substances
line1 =fgetl(f1);
y = strsplit(line1);
l =length(y);
ele =y(1);
elements = char(ele);
element=['H', 'C', 'O', 'N', 'Ar'];
atomic_weight = [1.008 12.011 15.999 14.007 39.948];
lowtemp = str2double(y(l-3));
midtemp = str2double(y(l-2));
hightemp = str2double(y(l-1));
%now save the values of line2 of the element
line2 =fgetl(f1);
c = strfind(line2,'E');
% First line high temperature co-efficients
c1 = str2double(line2(1:c(1)+3));
c2 = str2double(line2(c(1)+4:c(2)+3));
c3 = str2double(line2(c(2)+4:c(3)+3));
c4 = str2double(line2(c(3)+4:c(4)+3));
c5 = str2double(line2(c(4)+4:c(5)+3));
%second line temperature values
line3 =fgetl(f1);
d = strfind(line3,'E');
c6 = str2double(line3(1:d(1)+3));
c7 = str2double(line3(d(1)+4:d(2)+3));
c8 = str2double(line3(d(2)+4:d(3)+3));
c9 = str2double(line3(d(3)+4:d(4)+3));
c10 = str2double(line3(d(4)+4:d(5)+3));
%third line temperature values
line4 =fgetl(f1);
e = strfind(line4,'E');
c11 = str2double(line4(1:e(1)+3));
c12 = str2double(line4(e(1)+4:e(2)+3));
c13 = str2double(line4(e(2)+4:e(3)+3));
c14 = str2double(line4(e(3)+4:e(4)+3));
LT =linspace(lowtemp,midtemp,100);
HT =linspace(midtemp,hightemp,100);
% Specific heat calculation
Cp = s_h(c1,c2,c3,c4,c5,c8,c9,c10,c11,c12,T,R,temp_mid);
% Enthalpy calculation
H = Enthalpy(c1,c2,c3,c4,c5,c6,c8,c9,c10,c11,c12,c13,T,R,temp_mid);
% Entropy calculation
S = Entropy(c1,c2,c3,c4,c5,c6,c8,c9,c10,c11,c12,c13,T,R,temp_mid);
% Calculating molecular weight of species
Molecularweight(i) = Molecular_weight(elements);
%FOLDER CREATION
file =elements;
mkdir(file);
current = pwd;
cd(file)
%creating the plot for specific heat
figure(1) % creating figure using figure command and name figure1
plot(T, Cp, 'linewidth', 3, 'color', 'r'); % plot the graph between T and Cp
% hold on
% plot(LT, SH_low, 'linewidth', 3, 'color', 'g');
% hold off
xlabel('Temperature (K)'); % xlabel
ylabel('Specific heat (KJ)'); % ylabel
% legend('LT','HT');
title(sprintf('specific heat vs Temperature range of %s', elements)); % Title for the graph
grid on % Turning the grid on
fig = sprintf('specific heat of %s.png',elements); % to depict the filename with name of the substance.png
saveas(gca,fig); % save figure as
%creating the plot for enthalpy
figure(2) % creating figure using figure command and name figure2
plot(T, H, 'linewidth', 3, 'color', 'g'); % plot the graph between T and H
% hold on
% plot(LT, H_low, 'linewidth', 3, 'color', 'g');
% hold off
xlabel('Temperature (K)'); % xlabel
ylabel('Enthalpy'); % ylabel
% legend('LT','HT');
title(sprintf('Enthalpy vs Temperature range of %s', elements)); % Title for the graph
grid on % Turning the grid on
fig = sprintf('Enthalpy of %s.png',elements); % to depict the file name with name of the substance.png
saveas(gca,fig); % save figure
% hold on
% creating the plot for entropy
figure(3) % creating figure using figure command and name figure2
plot(T, S, 'linewidth', 3, 'color', 'y'); % plot the graph between T and H
% hold on
% plot(LT, S_low, 'linewidth', 3, 'color', 'g');
% hold off
xlabel('Temperature (K)'); % xlabel
ylabel('Entropy'); % ylabel
% legend('LT','HT');
title(sprintf('Entropy vs Temperature range of %s', elements)); % Title for the graph
grid on % Turning the grid on
fig = sprintf('Entropy of %s.png',elements); % to depict the file name with name of the substance.png
saveas(gca,fig); % save figure
cd(current); % change directory to current
% Molecular weight
M = Molecular_weight(elements)
f2 = fopen('molecular mass','w');
fprintf('%s n %s n %f n','species','Molecular mass',M)
fclose(f2);
endThe main code will calculate the specific heat enthalpy and entropy based on the temperature values for various substances from the THERMO.dat file.
ERRORS:
1. Unrecognized function or variable T
Solution: The temperature range variable is defined as temperature and changed to T
2.Error using S_h, Too many input arguments
Solution: input arguments in the main program is modified to represent same arguments as function
3. Error using enthalpy, Too many input arguments
Solution: input arguments in the main program is modified to represent same arguments as function
4. Unrecognized function or variable SH_high
solution:in plot command changed the name of thje variable to SH_high to H
5. Unrecognized variable CP in line 100
Solution:change CP to Cp in plot command
Results:
1. Specific heat vs Temperature of O2
2. Enthalpy vs temperature of O2
3. Entropy vs Temperature of O2
1. Specific heat vs Temperature of N2
2. Enthalpy vs temperature of N2
3. Entropy vs Temperature of N2
1. Specific heat vs Temperature of CO2
2. Enthalpy vs temperature of CO2
3. Entropy vs Temperature of CO2
Summary:
The project was done starting from extracting the data from the THERMO.dat file. We have developed the Cp, Enthalpy, and Entropy functions based on the equations given to us in the project description after selecting the temperature range from low to high with 1000 values and plotting the three functions vs temperature to get the results.
The error files and location of output files with the code file locations screenshots are attached below. The molecular weight screenshots from the command window are also attached below.
References:
1.Herbert B callen Thermodynamics an introduction to physical theories of equilibrium thermostatics
2.https://www.mathworks.com/help/matlab/ref/fopen.html
3.https://www.mathworks.com/help/matlab/ref/fgetl.html
4.https://www.mathworks.com/help/matlab/ref/mkdir.html
5.https://www.mathworks.com/help/matlab/ref/save.html
6.https://www.mathworks.com/help/matlab/ref/saveas.html
7.https://www.mathworks.com/matlabcentral/answers/353089-how-to-solve-the-fopen-error
8.https://www.mathworks.com/matlabcentral/answers/1803135-error-with-fgetl-why