Parsing NASA thermodynamic data
Objective: To Parse the thermodynamic data and use the data to make plots and generate information such as the molecular weight.
Problem Statement:
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.
To guide you we have broken down the problem into meaningful bits that you can work on individually.
1. Write a function that extracts the 14 coefficients 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 temperature
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.
Program:
clear all
close all
clc
f1 = fopen('THERMO.dat','r');
% Reading First Line
fgetl(f1);
% Scaning numerical values
a = fscanf(f1,'%f');
global_low_temp = a(1);
global_mid_temp = a(2);
global_high_temp = a(3);
% Skip next three lines
for i = 1:3
fgetl(f1);
end
% For each element the plots and data are generated through for loop
for j = 1:53
species = fgetl(f1);
species_name = species(1:6);
z=49;
local_low_temp = str2num(species(z:z+8))
local_high_temp = str2num(species(z+9:z+17))
local_mid_temp = str2num(species(z+18:z+26))
line_1 = fgetl(f1);
c = findstr(line_1,'E');
a1 = str2num(line_1(1:c(1)+3));
a2 = str2num(line_1(c(1)+4:c(2)+3));
a3 = str2num(line_1(c(2)+4:c(3)+3));
a4 = str2num(line_1(c(3)+4:c(4)+3));
a5 = str2num(line_1(c(4)+4:c(5)+3));
line_2 = fgetl(f1);
c_1 = findstr(line_2,'E');
a6 = str2num(line_2(1:c_1(1)+3));
a7 = str2num(line_2(c_1(1)+4:c_1(2)+3));
b1 = str2num(line_2(c_1(2)+4:c_1(3)+3));
b2 = str2num(line_2(c_1(3)+4:c_1(4)+3));
b3 = str2num(line_2(c_1(4)+4:c_1(5)+3));
line_3 = fgetl(f1);
c_2 = findstr(line_3,'E');
b4 = str2num(line_3(1:c_2(1)+3));
b5 = str2num(line_3(c_2(1)+4:c_2(2)+3));
b6 = str2num(line_3(c_2(2)+4:c_2(3)+3));
b7 = str2num(line_3(c_2(3)+4:c_2(4)+3));
% Creating data for storing the temperature between a range
tempr = linspace(local_low_temp,local_high_temp,50);
% Universal Gas Constant
r = 8.314;
% Molecular Weight Function call
molecular_weight = mol_weight(species_name);
disp('Molecular weight of ');
disp(species_name);
disp('is ');
disp(molecular_weight);
% Creating directory to store the plots
mkdir(['C:\Users\narenderkumar\Desktop\ju\Mechanical Course\MatLab\Challenge 9\',species_name]);
% Based on the value being greater than mid value or lower than that
% calculations are performed
for k = 1:length(tempr)
t = tempr(k);
if tempr(k)<=local_mid_temp
cp(k) = (b1 + b2.*t + b3.*(t.^2)+b4.*(t.^3)+b5.*(t.^4)).*r;
h(k) = (b1 + (b2.*t)./2 + (b3.*(t.^2))./3 + (b4.*(t.^3))./4 + (b5.*(t.^4))./5 + b6./t)*r.*t;
s(k) = (b1.*log(t) + b2.*t + (b3.*(t.^2))./2 + (b4.*(t.^3))./3 + (b5.*(t.^4))./4 + b7).*r;
else
cp(k) = (a1 + a2.*t + a3.*(t.^2)+a4.*(t.^3)+a5.*(t.^4)).*r;
h(k) = (a1 + (a2.*t)./2 + (a3.*(t.^2))./3 + (a4.*(t.^3))./4 + (a5.*(t.^4))./5 + a6./t)*r.*t;
s(k) = (a1.*log(t) + a2.*t + (a3.*(t.^2))./2 + (a4.*(t.^3))./3 + (a5.*(t.^4))./4 + a7).*r;
end
end
% Plotting the graph
figure(1)
plot(tempr,cp)
xlabel('Temperature')
ylabel('Specific Heat')
title(species_name)
figure(2)
plot(tempr,h)
xlabel('Temperature')
ylabel('Enthalpy')
title(species_name)
figure(3)
plot(tempr,s)
xlabel('Temperature')
ylabel('Entropy')
title(species_name)
% Saving files by changing the directory to the exisiting element
dir = (['C:\Users\narenderkumar\Desktop\ju\Mechanical Course\MatLab\Challenge 9\',species_name]);
cd(dir);
saveas(1,'Specific_Heat.jpg');
saveas(2,'Enthalpy.jpg');
saveas(3,'Entropy.jpg');
cd('C:\Users\narenderkumar\Desktop\ju\Mechanical Course\MatLab\Challenge 9')
end
fclose(f1)
Steps:
- First, we write all the clear statements to remove any previous trash data.
- Now, we use the fopen function to open the data file provided i.e 'THERMO.dat' and in fopen mention 'r' which means read.
- Then we obtain the first line by using the function fgetl(f1) where the f1 is already assigned with the opened file. In the next step, we obtain all the values for global low, mid, and high using the fscanf where only numeric data is extracted.
- Then we use the for loop to skip the next 3 lines.
- In the next step, we create a for loop where we generate plots and data for each molecule.
- Here we use the fgetl once again to get the corresponding line and store in species and then obtain the name of the species.
- Later the value for z is mentioned based on the character position from where we can obtain the values of the temperatures and then store them in new variables.
- We repeat the same process for the following lines and obtain high and low coefficients using the findstr and then convert the data to numeric form using str2num
- A dataset of 50 is values is created between the lowest and highest temperatures and stored. the value of 'r' is assigned next.
- Then the function for molecular weight is run and the values are displayed, the function is created in the different file whose steps would be mentioned ahead.
- A directory is created using the mkdir and the location where the folders are to generated are mentioned along with their names.
- Then we run a for loop where using the formula we generate the values for specific heat, entropy, and enthalpy.
- Later, we plot the graph and then save these plots to the directory of the species.
Function for Molecular weight
function m_w = mol_weight(speciesname)
% Primary Atoms
elements = ['H','C','O','N','A'];
% Data for atomic_weights
atomic_weight = [1.0079 12.0107 15.099 14.0067 39.948];
m_w =0;
% Loop through specimen and then each element to get required data and
% process further
for i = 1:length(speciesname)
for j =1:length(elements)
if strcmp(speciesname(i),elements(j))
m_w = m_w + atomic_weight(j);
y = j;
end
end
n = str2num(speciesname(i));
if n>1
m_w = m_w + atomic_weight(y)*(n-1);
end
end
end
Steps:
- In this function, we write the code to calculate the molecular weight for each species and then displaying them
- Make sure that the name of the file and function are identical else, it will lead to an error.
- The main elements are stored in a variable in the array form along with their atomic weights in another variable.
- Then a for loop is written to perform the calculations for the molecular weight containing these elements in them and then finally giving the output for it.
- The strcmp is used to compare the strings which are mentioned above.
- Don't run this file as the input arguments are not passed on here.
Plots:
The folder containing all the elements and the plots in them.
Errors:
To avoid this error run the program from the main file and not where the function exists. Here we are not passing any input arguments so the MatLab throws an error.
Results:
These are a few of the molecular weights which are generated after running the main file.
References:
https://www.mathworks.com/help/matlab/ref/fscanf.html