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.

Code:

clear all
close all
clc

% To read the file
f1=fopen('THERMO.dat','r');
fgetl(f1); %To get the 1st line

%To get the second line
Temp=str2num(fgetl(f1)); % To convert string to number
global_low_temp=Temp(1);
global_mid_temp=Temp(2);
global_high_temp=Temp(3);

% Skipping unwanted lines
for i=1:3
comments=fgetl(f1);
end

% For loop to read the values and plot molecular weight
for j=1:53
species=fgetl(f1);
S=strsplit(species);
species_name=(S{1});

P=length(S);
local_low=str2double(S{P-3});
local_high=str2double(S{P-2});
local_mid=str2double(S{P-1});

line_1=fgetl(f1);
A=findstr(line_1,'E');
a1=str2num(line_1(1:(A(1)+3)));
a2=str2num(line_1((A(1)+4):(A(2)+3)));
a3=str2num(line_1((A(2)+4):(A(3)+3)));
a4=str2num(line_1((A(3)+4):(A(4)+3)));
a5=str2num(line_1((A(4)+4):(A(5)+3)));

line2=fgetl(f1);
B=findstr(line2,'E');
a6=str2num(line2(1:(B(1)+3)));
a7=str2num(line2((B(1)+4):(B(2)+3)));
a8=str2num(line2((B(2)+4):(B(3)+3)));
a9=str2num(line2((B(3)+4):(B(4)+3)));
a10=str2num(line2((B(4)+4):(B(5)+3)));

line3=fgetl(f1);
c=findstr(line3,'E');
a11=str2num(line3(1:(c(1)+3)));
a12=str2num(line3((c(1)+4):(c(2)+3)));
a13=str2num(line3((c(2)+4):(c(3)+3)));
a14=str2num(line3((c(3)+4):(c(4)+3)));

% Creating temperature values temperature
temp=linspace(local_low,local_high,100);

% Value of the initial gas constant in J/mol*K
R=8.314;

% Calling specific heat function, enthalpy function and entropy function
cp=specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,R,temp,local_mid);
h=enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,R,temp,local_mid);
s=entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,R,local_mid);

figure(1)
plot(temp,cp)
xlabel('Temperature')
ylabel('Specific Heat')
title('Temperature vs Specific Heat')

figure(2)
plot(temp,h)
xlabel('Temperature')
ylabel('Enthalpy')
title('Temperature vs Enthalpy')

figure(3)
plot(temp,s)
xlabel('Temperature')
ylabel('Entropy')
title('Temperature vs Entropy')

mkdir(species_name)
cd(species_name)
saveas(1,'Specific Heat','jpg')
saveas(2,'Enthalpy','jpg')
saveas(3,'Entropy','jpg')
cd ..

% Calling the molrcular weight function
Mw=molecular_weight(species_name);

end

Function Codes:

1)Specific Heat function code

%function program for specific heat

function[cp]=specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,R,temp,local_mid)

if temp<=local_mid
cp =R*(a8+(a9*temp)+(a10*(temp).^2)+(a11*(temp).^3)+(a12*(temp).^4));

else
cp =R*(a1+(a2*temp)+(a3*(temp).^2)+(a4*(temp).^3)+(a5*(temp).^4));
end
end

2) Enthalpy function code:

% Function program for enthalpy

function[h]=enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,R,temp,local_mid)
if temp<=local_mid
h=(a8+((a9*temp)/2)+((a10*(temp).^2)/3)+((a11*(temp).^3)/4)+((a12*(temp).^4)/5)+(a13./temp)).*temp*R;
else
h=(a1+((a2*temp)/2)+((a3*(temp).^2)/3)+((a4*(temp).^3)/4)+((a5*(temp).^4)/5)+a6./temp).*temp*R;

end
end

3)Entropy function code

%function program for entropy

function[s]=entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,R,local_mid)
if temp<=local_mid
s=R*(a8*log(temp)+a9*temp+((a10*((temp).^2)/2))+((a11*((temp).^3)/3))+((a12*((temp).^4)/4))+a14);

else
s=R*(a1*log(temp)+a2*temp+((a3*((temp).^2)/2))+((a4*((temp).^3)/3))+((a5*((temp).^4)/4))+a7);

end
end

4) Molecular weight function

%function program for molecular weight

function[Mw]=molecular_weight(species_name)

%elements
element={'H' 'C' 'N' 'O' 'Ar'};

%Atomic weight
atomic_weight=[1.00784 12.0107 14.0067 15.999 39.948];

Mw=0;

% For loop to calculate molecular weight
for i=1:length(species_name)
for j=1:length(atomic_weight)
if strcmp(species_name(i),element(j));
Mw=Mw+atomic_weight(j);

u=j;
end
end

n=str2num(species_name(i));
if n>1
Mw=Mw+(atomic_weight(u)*(n-1));
end
end
fprintf('Molecular weight of %s = %f',species_name,Mw);

end

Code explaination:

As always we started with

clear all

close all

clc

After that to read th file from 'Thermo.dat" we use fopen function and to get the first line from that file we use fgetl command

% To read the file
f1=fopen('THERMO.dat','r');
fgetl(f1); %To get the 1st line

after this to read second line which is the global temp line we also use fgetl command with str2num command to convert the string into the command

%To get the second line
Temp=str2num(fgetl(f1)); % To convert string to number
global_low_temp=Temp(1);
global_mid_temp=Temp(2);
global_high_temp=Temp(3);

Then we skip the 3 unwanted lines

% Skipping unwanted lines
for i=1:3
comments=fgetl(f1);
end

Now we are using for loop to specify the species and to plot the molecular weight

for j=1:53
species=fgetl(f1);
S=strsplit(species);
species_name=(S{1});

after that we declare the local temperatures with low,mid and high temp.

P=length(S);
local_low=str2double(S{P-3});
local_high=str2double(S{P-2});
local_mid=str2double(S{P-1});

Now to get the data from above line1,line2 and line3 we use following steps

line_1=fgetl(f1);
A=findstr(line_1,'E');
a1=str2num(line_1(1:(A(1)+3)));
a2=str2num(line_1((A(1)+4):(A(2)+3)));
a3=str2num(line_1((A(2)+4):(A(3)+3)));
a4=str2num(line_1((A(3)+4):(A(4)+3)));
a5=str2num(line_1((A(4)+4):(A(5)+3)));

line2=fgetl(f1);
B=findstr(line2,'E');
a6=str2num(line2(1:(B(1)+3)));
a7=str2num(line2((B(1)+4):(B(2)+3)));
a8=str2num(line2((B(2)+4):(B(3)+3)));
a9=str2num(line2((B(3)+4):(B(4)+3)));
a10=str2num(line2((B(4)+4):(B(5)+3)));

line3=fgetl(f1);
c=findstr(line3,'E');
a11=str2num(line3(1:(c(1)+3)));
a12=str2num(line3((c(1)+4):(c(2)+3)));
a13=str2num(line3((c(2)+4):(c(3)+3)));
a14=str2num(line3((c(3)+4):(c(4)+3)));

Now we creating temp values

% Creating temperature values temperature
temp=linspace(local_low,local_high,100);

Then we are giving the initial gas constant value R

R=8.314;

After that Calling specific heat function, enthalpy function and entropy function

cp=specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,R,temp,local_mid);
h=enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,R,temp,local_mid);
s=entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,R,local_mid);

Then we are plotting for specific heat, enthalpy and entropy

figure(1)
plot(temp,cp)
xlabel('Temperature')
ylabel('Specific Heat')
title('Temperature vs Specific Heat')

figure(2)
plot(temp,h)
xlabel('Temperature')
ylabel('Enthalpy')
title('Temperature vs Enthalpy')

figure(3)
plot(temp,s)
xlabel('Temperature')
ylabel('Entropy')
title('Temperature vs Entropy')

After that we are creating the folders to save the plots

mkdir(species_name)
cd(species_name)
saveas(1,'Specific Heat','jpg')
saveas(2,'Enthalpy','jpg')
saveas(3,'Entropy','jpg')
cd ..

After that we calling the molecular weight function

Mw=molecular_weight(species_name);

At the end we put 'end' to end the for loop.

Plots;

O2:

N2:

CO2: