Project

Parsing NASA Thermodynamic Data using MATLAB

KKhaled Ali12 Dec 2019

MATLAB

Objective:

write a program that extracts the 14 coefficients from given THERMO.DAT file
Calculate the enthalpy, entropy and specific heat for all the species in the data file
Calculate the molecular weight of each species and display it in the command window.
Plot the Cp, Enthalpy and Entropy for the local temperature range (specific for each species).
To save the plots as images with appropriate names and dump them in separate folders for each species.
To save the Cp, Enthalpy and Entropy data in a text file
To make a program to interact with the user in order to automate the whole process where the user enter the species name and the program then retrieves the appropriate data (the plots) and display it in front of the user.

Input Equations:

`(CP)/R=a1+a2*(T)+a3*(T)^2+a4*(T)^3+a5*(T)^4`

`H/(RT)=a1+(a2*T)/2+(a3*(T)^2)/3+(a4*(T)^3)/4+(a5*(T)^4)/5+(a6)/T`

`S/R=a1*ln(T)+a2*T+(a3*(T)^2)/2+(a4*(T)^3)/3+(a5*(T)^4)/4+a7`

Molecular weight of a species is calculated by summing the atomic weights

Input:

`R = 8.3145 J/(K-mol)`

Atomic weights of the atoms are given below :

`H =1.00784`

`C = 12.0107`

`Ar = 39.948`

`N =14.0067`

`O =15.999`

Main code:

clear all
close all
clc

R=8.314;
f1=fopen('THERMO.DAT','r');
%preprocessing 
header=fgetl(f1);
overalglobaltemperatureranges=fgetl(f1);
a=strsplit(overalglobaltemperatureranges,' ');
low_global_temperature=str2num(a{2});
medium_global_temperature=str2num(a{3});
high_global_temperature=str2num(a{4});
for i=1:3;
    erase=fgetl(f1);
end
%parsing the species
for j=1:53;
 overal_local_temperature_ranges=fgetl(f1);
b=strsplit(overal_local_temperature_ranges,' ');
low_local_temperature=str2num(b{end-3});
medium_local_temperature=str2num(b{end-1});
high_local_temperature=str2num(b{end-2});
species_name=b{1};
temperature_range=linspace(low_local_temperature,high_local_temperature,200);
coef_line_1=fgetl(f1);
coef_line_2=fgetl(f1);
coef_line_3=fgetl(f1);
    % Calculating the Values of cofficients
x=findstr(coef_line_1,'E');
a1=str2num(coef_line_1(1:x(1)+3));
a2=str2num(coef_line_1(x(1)+4:x(2)+3));
a3=str2num(coef_line_1(x(2)+4:x(3)+3));
a4=str2num(coef_line_1(x(3)+4:x(4)+3));
a5=str2num(coef_line_1(x(4)+4:x(5)+3));

y=findstr(coef_line_2,'E');
a6=str2num(coef_line_2(1:y(1)+3));
a7=str2num(coef_line_2(y(1)+4:y(2)+3));
a8=str2num(coef_line_2(y(2)+4:y(3)+3));
a9=str2num(coef_line_2(y(3)+4:y(4)+3));
a10=str2num(coef_line_2(y(4)+4:y(5)+3));

z=findstr(coef_line_3,'E');
a11=str2num(coef_line_3(1:z(1)+3));
a12=str2num(coef_line_3(z(1)+4:z(2)+3));
a13=str2num(coef_line_3(z(2)+4:z(3)+3));
a14=str2num(coef_line_3(z(3)+4:z(4)+3));


    % Calculating the Values of Cp, Enthalpy and Entropy
[cp h s] = thermodynamic_properties(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,temperature_range,medium_global_temperature);
     % Calculating the Molecular Weight Of Species
[molecular_weight(j)]=molecular_mass(species_name);

 %plotting
set(0,'DefaultFigureVisible','off')
m(1)=figure;
plot(temperature_range,cp,'linewidth',3,'color','g');
xlabel('Temperature (K)');
ylabel('Specific Heat (J/kg-K)');
title(['Specific Heat vs Temperature ',species_name]);
 
 
 m(2)=figure ; 
 plot(temperature_range,h,'linewidth',3,'color','r');
 xlabel('Temperature (K)');
 ylabel('Enthalpy (J/kg)');
 title(['Enthalpy vs Temperature ',species_name]);
 
 m(3)=figure ;
 plot(temperature_range,s,'linewidth',3,'color','y');
 xlabel('Temperature (K)');
 ylabel('Entropy (J/kg-K)');
 title(['Entropy vs Temperature ',species_name]);
 
 % Printing Molecular Weight Data in Command Window
    fprintf ('Species Name : %6s     ',species_name)
    fprintf ('Molecular Weight : %15s \n',molecular_weight(j))
   
    cd('D:\species');
    mkdir(species_name);
    cd(species_name);
    % Saving as .jpg for windows explorer
    saveas(m(1),'specific_heat_vs_Temp.jpg');
    saveas(m(2),'enthalpy_vs_Temp.jpg');
    saveas(m(3),'entropy_vs_Temp.jpg');
    % Saving as .fig for MATLAB
    saveas(m(1),'specific_heat_vs_Temp.fig');
    saveas(m(2),'enthalpy_vs_Temp.fig');
    saveas(m(3),'entropy_vs_Temp.fig');

    
    %saving data to file
    f2 = fopen('data_parsing.txt','w');
    fprintf (f2, '%15s : %15s \n','Species Name',species_name);
    fprintf (f2, '%15s : %15.3s \n','MolecularWeight',molecular_weight(j));
    fprintf (f2, 'For the Local Temeprature Range : %3.1fK to %3.1fK \n\n',low_local_temperature, high_local_temperature);
    fprintf (f2, '%6s %6s %16s %16s \n','Temperature(K)','Specific Heat(J/K)','Enthalpy(J/kg)','Entropy(J/K)');
    fprintf (f2, '%11.2f %17.2f %19.2f %20.2f \n',temperature_range,cp,h,s);
    
    
    cd ('D:\Matlab\bin');
end


fprintf('\n The Input File has been stored \n\n')

% Getting input From User To Display Plots 
    Input = input('Do you require Data for Specific Species? Y/N : \n ','s')

if strcmpi(Input,'Y')
    Name = input ('Provide Name of Species : ','s')
    cd('D:\species')
    cd (Name)
    openfig('specific_heat_vs_Temp','visible');
    openfig('enthalpy_vs_Temp','visible');
    openfig('entropy_vs_Temp','visible');
    fprintf ('The Plots Corresponding to the Species have been Opened \n')
    fprintf ('\n The Program has ended \n')
    
elseif strcmpi(Input,'N')
    fprintf ('\n The Program has ended \n')
    
else
    fprintf ('Invalid Input \n')
    fprintf ('\n The Program has ended \n')
end

Main code explaination:

First input R value is declared
Then the file Thermo.Dat is read
Then Fgetl command is used to read the next line and then it is stored in header
Then the next line is read and stored in overall global temperatures ranges
Then strsplit command is used to split the string using the space as a point of split
Then str2num command is used to convert the string into a number and the second and the third and the fourth element of the converted string is stored as global low ,mid, and high temperatures respectively and the cell brackets are used{}
Then for loop is used with fgetl to read and store the next 3 lines in erase.
Then a loop for parsing elements are made and its counter is set to 53 because it is the number of elements
Then fgetl is used to store the next line and store it in overall local temperature ranges
Then strsplit is used to split the overall local temperature string using the space as a point of split and store it in 'b'
Then this string is converted into number with its first element before the end and second element before the end and third element before the end stored in medium, high, and low local temperatures respectively, and the cell brackets are used{}.
Then the first element in cell b is stored in species name as character
Then linspace command was used in order to specify the temperature range on the x axis and 200 values was selected as the number of included values including the bounds
Then the three lines of coefficients is read by fgetl command and stored
Then findstr command is used in order to find the location of the letter 'E' and it is stored in the value x
Then all the places before this location which is stored in x are selected and three places after as this is the trend in all coefficients in order to capture all the coefficient value
Then the next value is specified using the first x location in addition to 4 places to 3 places after the second location stored in x. This is done to all coefficients in order to capture all the coefficient values.
Then this is repeated for all coefficient and the 'E' locations for second and third lines are stored in y and z respectively.
Then function thermodynamic_properties is called and the output is specified and its input specified
Then molecular_mass function is called and its output specified in an array using ’j’ as a counter
Then figure visibility is set to off
Then figure is stored in m
Then all the plots are plotted using plot command and linewidth, color of the plot , axle labels and graph title are specified.
Then fprintf command is used to print molecular weight data in the command window and data type string is specified and the space is specified and species_name and molecular weight array are selected as inputs.
Then MATLAB is ordered to change directory to new directory species
Then inside it, it will make a directory for each species
Then it will go inside (change)to this directory
Then it will save as all the plots which were plotted as jpg and fig to open in windows explorer and matlab
Then using fopen command ,MATLAB will save (write) test data to data_parsing file
Then all the data to be written are selected and spaces and its general formatting is specified
Then MATLAB returns to its original location to complete the program
Then the loop is ended
The using fprintf,'the input file has been stored' is printed in command screen
Then user program starts
Using the Input command, the user is asked if it needs any species plots
Then strcmpi command is used to compare user answer to input irrespective of it being lower-case or upper-case letter
If user said 'y' ,the program will ask to provide its name which will be stored in name
Then matlab changes directory to main folder
And then changes directory to name which is the same as species name which store all the plots for any species
Then all the required plots are opened and made visible
Then using fprinf ,a message that the corresponding plots are opened is printed in command screen
Then using fprintf the message 'the program has ended' is written in command window
Then If user said 'N' , using fprintf the message 'the program has ended' is written in command window
Then If user entered wrong input, 'invalid input message' will appear in command screen using fprintf
Then using fprintf the message 'the program has ended' is written in command window

Thermodynamic properties function code:

function[cp,h,s]=thermodynamic_properties(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,temperature_range,medium_global_temperature)
for k=1:length(temperature_range);
if temperature_range(k)>medium_global_temperature
    cp(k)=R.*(a1+(a2.*temperature_range(k))+(a3.*(temperature_range(k)).^2)+(a4.*(temperature_range(k)).^3)+(a5.*(temperature_range(k).^4)));
    h(k)=R.*temperature_range(k).*(a1+(a2.*temperature_range(k))./(2)+(a3.*(temperature_range(k)).^2)./(3)+(a4.*(temperature_range(k)).^3)./(4)+(a5.*(temperature_range(k)).^4)./(5)+(a6)./(temperature_range(k)));
    s(k)=R.*((a1.*log(temperature_range(k)))+(a2.*temperature_range(k))+(a3.*(temperature_range(k)).^2)./(2)+(a4.*(temperature_range(k)).^3)./(3)+(a5.*(temperature_range(k)).^4)./(4)+(a7));
    
    else
    cp(k)=R.*(a8+(a9.*(temperature_range(k))+(a10.*(temperature_range(k)).^2)+(a11.*(temperature_range(k)).^3)+(a12.*(temperature_range(k)).^4)));
    
    h(k)=R.*temperature_range(k).*(a8+(a9.*temperature_range(k))./(2)+(a10.*(temperature_range(k)).^2)./(3)+(a11.*(temperature_range(k)).^3)./(4)+(a12.*(temperature_range(k)).^4)./(5)+(a13)./(temperature_range(k)));
    s(k)=R.*((a8.*log(temperature_range(k)))+(a9.*temperature_range(k))+(a10.*(temperature_range(k)).^2)./(2)+(a11.*(temperature_range(k)).^3)./(3)+(a12.*(temperature_range(k)).^4)./(4)+a14);
end
end

Thermodynamic properties function code Explaination:

First the function output and input are specified
Then using k loop which starts at 1 and ends at the length of temperature range input, the equations are calculated for each temperature value provided from the temperature range array
Then using 'if' ,a specific temperature is picked from the array and it is compared to the medium global temperature in order to specify if the lower or high temperature coefficient that will be used
Then all the outputs are put in an array using counter k in order to store 200 value for each species to be able to then plot it

Molecular weight Function code:

function [molecular_weight] = molecular_mass(species_name)
atoms=['O','C','H','A','N'];
atomic_weight=[15.999,12.0107,1.00784,39.948,14.0067];
molecular_weight=0;
x=1;
for i=1:length(species_name);
for j=1:length(atoms);
        
   if strcmp(species_name(i),atoms(j))==1 ;
      molecular_weight=atomic_weight(j)+molecular_weight;
      atomicweight=atomic_weight(j);
   end
end
   t=isstrprop(species_name(i),'digit');
   c=str2double(regexp(species_name,'\d*','Match'));     
   
 if t==1
    molecular_weight=atomicweight*(c(x)-1)+molecular_weight;
    x=x+1;
   
 end 
 
end
end

Molecular weight Function code Explaination:

First the function output is declared and function input which is a character is selected
Then array of atom characters is declared
Then array of atomic weight is declared with the same arrangement as the atom characters array
Then initial molecular weight is declared
Then initial x value is declared
The 'i' counter is set by the input species_name length
Then 'j' counter is set by array of atom character length
So, each input is compared to all the atom array using strcmp command
If this comparison is positive the program calculates the atomic weight by adding the old molecular weight to the atomic weight of the current species which is specified by the counter j because the element in atoms array will have the same j value as its atomic weight array because they are arranged in the same sequence
Then isstrprop command is used to test if this string is digit or not and it outputs a logical value either 0 or 1 and it is stored in 't'
Then using regexp command to extract all the digits from the input species_name and then they are converted from a cell array of character vectors to a numeric array using str2double command.
Then the program performs another test using if condition when t equals one which means that this string is a digit
Then the molecular weight will be equal the old molecular weight in addition to the c value which contains the digit value (picked by counter x) minus one because already the program has calculated the molecular weight once for the element before it had moved over to the digit
Then x which is used as a counter is increased by one in order when the next digit comes ,the program will select the right value in the the c array and won’t select the all the elements together.

Errors during programing:

This error happening because of using " sign instead of ' sign.It was then corrected.

This error happened because of the array that was produced contain more than one element in c variable,so it was corrected by implementing counter x and then using it to pick certain element from the array

This error happened because of wrong equation brackets it was then corrected

Also there was the problem of plots popping up while the program as running which was corrected using set(0,'DefaultFigureVisible','off') command and assigning each figure to variable.

OUTPUT:

O2 PLOTS

CO2 PLOTS

N2 PLOTS

MOLECULAR WEIGHT PRINT IN COMMAND WINDOW

TEXT FILE OUTPUT

FOLDER DIRECTORY

ALL OUTPUT FOLDERS LINK:

All output folders can be found at this link:

https://drive.google.com/drive/folders/1mnOgNaqogiySR0p1-6cVNesSx3yWFi8m?usp=sharing

References:

https://www.tutorialspoint.com/matlab/matlab_strings.htm