Challenge

Project 1 - Parsing NASA thermodynamic data

SShreyas Dabhade28 Sept 2022

Aim:

To write a Matlab program and calculate the thermodynamic properties of various gas species by parsing the NASA thermodynamic data file.

Objective:

To write a function that extracts the 14 coefficients and calculates the enthalpy, entropy, and specific heats for all the species in the data file.
Calculate the molecular weight of each species.
Plot the Specific heat, Enthalpy, and Entropy for the local temperature range (low temperature: high temperature) specific for each species.
To save the plots as images with appropriate names and dump them in separate folders for each species with a suitable name.

Theory:

File Parsing:

It is the method of splitting the file or other input into a piece of data that can be easily stored or manipulated. By doing this we are going to split a string into parts and then recognize the parts to convert it into something simpler than a string.

Here in this project, we are going to split the NASA thermodynamic data which is having 53 gas species values in 213 lines. To perform the calculation for thermodynamic properties, we need some sort of co-efficient values and temperature ranges which is available in the NASA file. So, we are going to parse the values which are all needed for our calculation from that file.

Specific Heat:

Specific heat is the amount of heat energy required to raise the temperature of the body per unit of mass. Specific heat is also known as specific heat capacity or mass-specific heat. In SI units, specific heat is the amount of heat in joules required to raise 1 gram of a substance by 1 Kelvin.

Enthalpy:

When a substance changes at constant pressure, enthalpy tells how much heat and work was added or removed from the substance. Enthalpy is similar to energy but not the same. When a substance grows or shrinks, energy is used up or released.

Entropy:

The entropy of an object is a measure of the amount of energy that is unavailable to do work. Entropy is also a measure of the number of possible arrangements of atoms a system can have. In this sense, entropy is the measure of uncertainty or randomness.

Gas Constant:

The gas constant is a physical constant denoted by R and is expressed in terms of units of energy per temperature increment per mole. The gas constant value is equivalent to the Boltzmann constant but expressed as the pressure-volume product instead of energy per increment of temperature per particle.

Molecular Weight:

Molecular weight is a measure of the sum of the atomic weight values of the atoms in a molecule. Molecular weight is used in chemistry to determine stoichiometry in chemical reactions and equations.

Formula:

To calculate the thermodynamic properties:

Here,

Cp = Specific Heat (J/K kg)

H = Enthalpy (J)

S = Entropy (J/K)

R = Universal gas constant = 8.314 (J/mol.K)

T = Temperature (K)

Sample NASA data:

1. Code FOR SPECIFIC HEAT ;

function [Cp]=specific_heat1(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,local_mid_temp)

% Calculation of specific heat for higher and lower temperature

if (T>=local_mid_temp)

Cp=R*(a1+(a2*T)+(a3*T).^2)+(a4*(T).^3)+(a5*(T).^4);

else

Cp=R*(a8+(a9*T)+(a10*T).^2)+(a11*(T).^3)+(a12*(T).^4);

end

Explanation of the above code: -

%Funtion for the Specific Heat - comments to make a program user friendly
function [Cp] = Specific_heat(x1,x2,x3,x4,x5,x8,x9,x10,x11,x12,temp,gb_m,r) - function is used to make the code easy to implement ,it is a group of statements that together perform a task
%Checking for the temperature -comments
%The FIRST 7 coefficients are HIGH-temperature coefficients and the SECOND 7 coefficients are LOW-temperature coefficients
if(temp<=gb_m) - comparing matrix values of the temperature generated from the linspace command with the global medium temperature
Cp = r.*(x8 +(x9.*temp) +(x10*(temp).^2) + (x11*(temp).^3) + (x12*(temp).^4)); - calculation of the Specific heat wrt to the equations
else - else statement to perform the task if the if statements fail
Cp = r.*(x1 +(x2.*temp) +(x3*(temp).^2) + (x4*(temp).^3) + (x5*(temp).^4)); - calculation of the Specific heat wrt to the equations
end

2. Entropy

function [H]=Enthalpy2(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,R,local_mid_temp)

% Calculation of Enthalpy for higher and lower temperature

if (T>=local_mid_temp)

H=R*((a1*T)+((a2*(T).^2)/2)+((a3*(T).^3)/3)+((a4*(T).^4)/4)+((a5*(T).^5)/5)+a6);

else

H=R*((a8*T)+((a9*(T).^2)/2)+((a10*(T).^3)/3)+((a11*(T).^4)/4)+((a12*(T).^5)/5)+a13);

end

Explanation of the above code: -

%Funtion for the entropy - comments to make a program user friendly
function [S] = entropy_S(x1,x2,x3,x4,x5,x7,x8,x9,x10,x11,x12,x14,temp,gb_m,r) -

the function is used to make the code easy to implement,it is a group of statements that together perform a task
%Checking for the temperature -comments
%The FIRST 7 coefficients are HIGH-temperature coefficients and the SECOND 7 coefficients are LOW-temperature coefficients - comments
if(temp<=gb_m) - comparing matrix values of the temperature generated from the linspace command with the global medium temperature
S = r.*(x8.*log(temp) +(x9.*temp) +(x10*((temp).^2)/2) + (x11*((temp).^3)/3) + (x12*((temp).^4)/4) + (x14)); - calculation of the entropy wrt to the equations
else - else statement to perform the task if the if statements fail
S = r.*(x1.*log(temp) +(x2.*temp) +(x3*((temp).^2)/2) + (x4*((temp).^3)/3) + (x5*((temp).^4)/4) + (x7)); - calculation of the entropy wrt to the equations
end - to end or terminate the function

3. Enthalpy:

function [S]=Entropy3(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,local_mid_temp)

% Calculation of Entropy of higher and lower temperature

if(T>=local_mid_temp)

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

else

S=R*(a8*log(T)+(a9*T)+((a10*(T).^2)/2)+((a11*(T).^3)/3)+((a12*(T).^4)/4)+a14);

end

Explanation of the above code: -

%Funtion for the enthalpy - comments to make a program user friendly
function [H] = enthalpy_H(x1,x2,x3,x4,x5,x6,x8,x9,x10,x11,x12,x13,temp,gb_m,r) - function is used to make the code easy to implement ,it is a group of statements that together perform a task
%Checking for the temperature -comments
%The FIRST 7 coefficients are HIGH-temperature coefficients and the SECOND 7 coefficients are LOW-temperature coefficients - comments
if(temp<=gb_m) - comparing matrix values of the temperature generated from the linspace command with the global medium temperature
H = r*temp.*(x8 + ((x9*temp)/2) + ((x10*(temp).^2)/3) + ((x11*(temp).^3)/4) + ((x12*(temp).^4)/5) + ((x13./temp))); - calculation of the enthalpy wrt to the equations
else - else statement to perform the task if the if statements fail
H = r*temp.*(x1 + ((x2*temp)/2) + ((x3*(temp).^2)/3) + ((x4*(temp).^3)/4) + ((x5*(temp).^4)/5) + ((x6./temp))); - calculation of the enthalpy wrt to the equations
end - to terminate the function

4. Molecular weight function:

function Mol_wt=Molecular_weight(species_name);

Atomic_names=['H' 'C' 'O' 'N' 'A' 'S'];

Atomic_weight=[1 12 16 14 40 32];

Mol_wt=0;

for i=1:length(species_name)

for j=1:length(Atomic_names)

if strcmpi(species_name(i),Atomic_names(j))

Mol_wt=Mol_wt+Atomic_weight(j);

position=j;

end

n=str2double(species_name(i));

if n>1

Mol_wt=Mol_wt+Atomic_weight(position)*(n-1);

end

Explanation of the above code: -

%function calculate molecular mass - comments to make a program user friendly
function mw = mol_wt(sp_name)
%storing value - comments to make a program user friendly
compound_name = (sp_name)
key_at = ['C','H','O','N','A']; - creating an array of elements
%storing the atomic mss value in the array - comments to make a program user friendly
value_at = [12.0107 1.00794 15.994 14.0067 39.948]; - storing the atomic mass value in the array

mw =0 ; -assigning molecular mass as 0
%creating for loop to perform the task - comments to make a program user friendly
for i = 1:length(compound_name) - creating for loop to execute till the elements length
for j = 1:length(key_at) - creating for loop to execute till the key_at array length
%comparing the string value - comments to make a program user friendly
if strcmp(compound_name(i),key_at(j)) - comparing the string value of the species name and the elements in the array key_at
%if the above comparison is the true value of molecular mass is - comments
%calculated- comments
mw =mw + value_at(j); - molecular mass calculation
%storing the position for further use- comments
pos = j; - storing the position
end - to terminate the if statement
end - to terminate the inner for j loop
k = str2num(compound_name(i)); - converting the string to number
%comparing if the value of k is greater than executed the further statements- comments
if(k>1)
%molecular mass calculation mw = mw + value_at(pos)*(k-1);- comments
end - to terminate the if statement
end- to terminate the outer for i loop
%command to print the molecular mass of the molecule- comments
fprintf('Molecular Weight of the species %s is %f',compound_name ,mw) - to print the species name and value
end - to terminate the function

Main Code :

% To write a matlab program for parsing nasa thermodynamic file

clear all

close all

clc

f1=fopen('THERMO.dat','r'); % Opening a file thermo.dat stored in memory

file=fgetl(f1)

local_temp=fgetl(f1)

% Gas constant

R=8.314;

% Converting string cells into numbers

local_temp_range=strsplit(local_temp,' ');

local_low_temp=str2double(local_temp_range{2})

local_mid_temp=str2double(local_temp_range{3})

local_high_temp=str2double(local_temp_range{4})

% Temperature range

T=linspace(local_low_temp,local_high_temp,1000);

% Skipping the lines

for i =1:3;

a=fgetl(f1);

end

for k = 1:53;

line1=fgetl(f1);

A=strsplit(line1,' ');

species_name=A{1};

line_1=fgetl(f1);

b=strfind(line_1,'E');

% Extracting 7 high temp coefficients and converting string to number

a1=line_1(1:b(1)+3);

a1=str2double(a1);

a2=line_1(b(1)+4:b(2)+3);

a2=str2double(a2);

a3=line_1(b(2)+4:b(3)+3);

a3=str2double(a3);

a4=line_1(b(3)+4:b(4)+3);

a4=str2double(a4);

a5=line_1(b(4)+4:b(5)+3);

a5=str2double(a5);

line_2= fgetl(f1);

c=strfind(line_2,'E');

a6=line_2(1:c(1)+3);

a6=str2double(a6);

a7=line_2(c(1)+4:c(2)+3);

a7=str2double(a7);

% Calculating 7 lower order coefficients

a8=line_2(c(2)+4:c(3)+3);

a8=str2double(a8);

a9=line_2(c(3)+4:c(4)+3);

a9=str2double(a9);

a10=line_2(c(4)+4:c(5)+3);

a10=str2double(a10);

line_3=fgetl(f1);

d=strfind(line_3,'E');

a11=line_3(1:d(1)+3);

a11=str2double(a11);

a12=line_3(d(1)+4:d(2)+3);

a12=str2double(a12);

a13=line_3(d(2)+4:d(3)+3);

a13=str2double(a13);

a14=line_3(d(3)+4:d(4)+3);

a14=str2double(a14);

% Calculation of specific heat

Cp=specific_heat1(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,local_mid_temp);

% Calculation of Enthalpy

H=Enthalpy2(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,R,local_mid_temp);

% Calculation of Entropy

S=Entropy3(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,local_mid_temp);

% Plotting T vs Cp

figure(1)

plot(T,Cp,'LineWidth',2,'color','g')

xlabel('Temperature[K]')

ylabel('Specific heat')

grid on

title('temperature vs specific heat')

% Plotting for T vs Enthalpy

figure(2)

plot(T,H,'LineWidth',2,'color','b')

xlabel('Temperature[K]')

ylabel('Enthalpy')

grid on

title('temperature vs Enthalpy')

% Plotting for T vs Entropy

figure(3)

plot(T,S,'LineWidth',2,'color','g')

xlabel('Temperature[K]')

ylabel('Entropy')

grid on

title('temperature vs Entropy')

% Molecular weight

Mol_wt=Molecular_weight(species_name);

k=sprintf('%s Mol_wt=%f',species_name,Mol_wt)

disp(['Molecular weight of: ',species_name])

disp(Mol_wt)

% Creating a directory

mkdir(['C:\Users\Admin\Desktop\matlab\Species1\',species_name])

cd(['C:\Users\Admin\Desktop\matlab\Species1\',species_name])

saveas(1,'cpvsT','jpg')

saveas(2,'HvsT','jpg')

saveas(3,'SvsT','jpg')

end

Explanation of the above CODE: -

%File Parsing for a file - comments to make a program user friendly
clc - to clear the command window
clear all - to clear all the workspace variables
close all - to close all the output window
f1 = fopen('THERMO.dat','r'); - command to open a file in read mode
first_line = fgetl(f1); - to read file line by line
%R is the universal gas constant J/(mol-K)
R = 8.314; - constant
temp = fgetl(f1); - to save the line in the temp variable

%Reading text from headerfile - comments to make a program user friendly
%Reading the global temperature from the next line
A = strsplit(temp,' '); - to split the line saved in temp variable when space(' ') is found

%Cell array is generated 1x4 cell array
glob_lowtemp = str2num(A{2}); - assigning the global lower temperature
glob_midtemp = str2num(A{3});- assigning the global middle temperature
glob_hightemp =str2num(A{4});- assigning the global high temperature

%Temperature range
Temperature = linspace(glob_lowtemp,glob_hightemp,1000); - command to generate the 1000 values between the range of temperature
%skipping the three comment lines
for j =1:3; - for skipping the comments line by running loop
k =fgetl(f1); - to read the comments line
end - end of the for loop
%running for loop fopr 53 times as there are 53 species
for i = 1:53; - for loop to parse the file and create plot for 53 molecules in the file

%Parsing of the species block -comments to make a program user friendly
tline = fgetl(f1); - reading the first line for a molecule
S =strsplit(tline,' '); - splitting the string when the space is found by the command
Species_name =S{1}; - extracting the name of the molecule

%Extracting the coefficients from first line
line2 =fgetl(f1); - reading the first line of the coefficients
x =findstr(line2,'E'); - finding the string using the command

%Finding the 5 coefficient
x1 = line2(1:x(1)+3); - Extracting the 1st coefficients from the line
x1 = str2num(x1); - converting the string coefficients extracted above to number
x2 = line2(x(1)+4:x(2)+3); - Extracting the 2nd coefficients from the line
x2 = str2num(x2); - converting the string coefficients extracted above to number
x3 = line2(x(2)+4:x(3)+3); - Extracting the 3rd coefficients from the line
x3 = str2num(x3); - converting the string coefficients extracted above to number
x4 = line2(x(3)+4:x(4)+3); - Extracting the 4th coefficients from the line
x4 = str2num(x4); - converting the string coefficients extracted above to number
x5 = line2(x(4)+4:x(5)+3); - Extracting the 5th coefficients from the line
x5 = str2num(x5); - converting the string coefficients extracted above to number
%line3 5coefficients
line3 = fgetl(f1); - reading the line of the coefficients
x = findstr(line3,'E');- finding the string using the command
x6 = line3(1:x(1)+3); -Extracting the 6th coefficients from the line
x6 = str2num(x6); - converting the string coefficients extracted above to number
x7 = line3(x(1)+4:x(2)+3);-Extracting the 7th coefficients from the line
x7 = str2num(x7); - converting the string coefficients extracted above to number
x8 = line3(x(2)+4:x(3)+3);-Extracting the 8th coefficients from the line
x8 = str2num(x8); - converting the string coefficients extracted above to number
x9 = line3(x(3)+4:x(4)+3);-Extracting the 9th coefficients from the line
x9 = str2num(x9); - converting the string coefficients extracted above to number
x10 = line3(x(4)+4:x(5)+3);-Extracting the 10th coefficients from the line
x10 = str2num(x10); - converting the string coefficients extracted above to number

%Line4 4coefficients
line4=fgetl(f1); - reading the line of the coefficients
x = findstr(line4,'E');- finding the string using the command
x11 = line4(1:x(1)+3);-Extracting the 11th coefficients from the line
x11 = str2num(x11); - converting the string coefficients extracted above to number
x12 = line4(x(1)+4:x(2)+3);-Extracting the 12th coefficients from the line
x12 = str2num(x12); - converting the string coefficients extracted above to number
x13 = line4(x(2)+4:x(3)+3);-Extracting the 13th coefficients from the line
x13 = str2num(x13); - converting the string coefficients extracted above to number
x14 = line4(x(3)+4:x(4)+3);-Extracting the 14th coefficients from the line
x14 = str2num(x14); - converting the string coefficients extracted above to number

% Entropy Calculation
S = entropy_S(x1,x2,x3,x4,x5,x7,x8,x9,x10,x11,x12,x14,Temperature,glob_midtemp,R);

%Cp calculation means Specific Heat calculation
Cp = Specific_heat(x1,x2,x3,x4,x5,x8,x9,x10,x11,x12,Temperature,glob_midtemp,R);

%Enthlpy Calculation is H
H = enthalpy_H(x1,x2,x3,x4,x5,x6,x8,x9,x10,x11,x12,x13,Temperature,glob_midtemp,R);

%Molecular Weight calculation of different elements
mw(i) = mol_wt(Species_name);

%Reading and writinh of the molecular wieght in the list of molecular
%wieght

% Plotting on the graph for species fo the temperature and specific
% heat

%Creating a folder for saving the species with its properties
directory_curr = pwd; - pwd identify the current directory
mkdir(Species_name); - command used for making a new directory
cd(Species_name); - cd is for changing the directory

%Plotting the graph in different plotting window
%plot 1 for the temperature and entropy
%plot 2 for the temperature and enthalpy
%plot 3 for the temperature and specific heat
%plot 1 Enthalphy and Temperaure
figure(1) - to open the first figure window
plot(Temperature,H,'linewidth',3,'color','b'); -plot between temperature and enthalpy
xlabel('Temperature in [K]'); - for labelling the x axis
ylabel('Enthalphy [KJ/mol]'); - for labelling the y axis
grid on - to make the grid on
title(sprintf('Variation in Enthalpy wrt Temperature %s',Species_name)); - to give the title to the plot
saveas (figure(1),'Enthalhy.jpg'); - to save the plot as jpg file in the folder made

%plot 2 Entropy and Temperature
figure(2) - to open the second figure window
plot(Temperature,S,'linewidth',3,'color','r'); -plot between temperature and entropy
xlabel('Temperature in [K]'); - labeling of the xaxis
ylabel('Entropy [KJ/mol-k]'); - labelling the y axis
grid on - to make the grid on
title(sprintf('Variation in Entropy wrt Temperature %s',Species_name));- to give the title to the plot
saveas (figure(2),'Entropy.jpg'); - to save the plot as jpg file in the folder made

%plot 3 Specific Heat and Temperature
figure(3) - - to open the third figure window
plot(Temperature,Cp,'linewidth',3,'color','g');plot between temperature and specific heat
xlabel('Temperature in [K]'); - labeling of the xaxis
ylabel('Specific Heat [KJ]');- labelling the y axis
grid on- to make the grid on
title(sprintf('Variation in Specific Heat wrt Temperature %s',Species_name));- to give the title to the plot
saveas (figure(3),'Specific Heat.jpg'); - to save the plot as jpg file in the folder made
cd(directory_curr); - to change the directory
end to end the for loop

Syntax Used during Coding

Fopen :	Open a file or obtain information about open files
‘r’ :	Open the file for reading
fgetl	Read line from file
strsplit	Split string or character vector at specified delimiter
str2num	Convert character array or string to numeric array

strfind	Find one string within another
linspace	Generate linearly spaced vectors
fprintf	Write formatted data to file
pwd	Identify current directory
mkdir	Make new directory
cd	Change working directory

Output / Result :

Attachment:

Screenshot of all the saved folders :

Results:

Specific Heat, Entropy, and Enthalpy Plots for O2, N2, and Co2 species :

1. Plots for O2

a. Variation of Specific heat with the temperature of O2

b. Variation of Entropy with the temperature of O2

c. Variation of Enthalpy with the temperature of O2

2. Plots for N2 :

a. Variation of Specific heat with the temperature of N2

b. Variation of Entropy with a temperature of N2

c. Variation of Enthalpy with the temperature of N2

3. Plots for Co2 :

a. Variation of Specific heat with the temperature of Co2

b. Variation of Entropy with the temperature of CO2

c. Variation of Enthalpy with the temperature of CO2

Conclusion :

Using a Matlab program, we have calculated the various coefficients required to evaluate various thermodynamic properties of the different species which are extracted from the NASA Thermodynamic data file.
The specific heat, Enthalpy, and Entropy are calculated for all the species and plots are plotted against the local temperature range and molecular weight for all the species is calculated.

Project 1 - Parsing NASA thermodynamic data

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.