MATLAB Project - Parsing NASA thermodynamic data

Problem Statement:-

1. Write a function that extracts the 14 co-effients and calculates the enthalpy, entropy and specifc heats for all the species in the data file.

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) specifc 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.

Theory :-

Problem Explanation :-

NASA came up with idea of polynomial which can be used to evaluate thermodynamic property such as `C_p , H and S`. With this idea NASA documented the co-efficients that are required to evaluate these polynomial for varies species.

In this problem we use such document which contain polynomial coefficient for 53 species. Screenshot of the file data is shown below.

This is the screeenshot of file. For better understanding some sections are highlighted with different color types.

File contains :-

1. First line in file(in sky blue) contain the name of file.
2. Second line Which is highlighted in Yello color contain Global temperature range. As first is global low temperature, second is global middle temperature and third is global ligh temperature.
3. Next three line which are highlighted in Sky blue color contain comment about file.
4. Next Four line are data of one species and consecutive next four line are data of other species and so on for other lines.
5. This four lines contains data of polynomial. In the right hand side of this lines, a number is given (which is highlighted in pink color), they define line number for particular species.
6. In species four lines, line 1 mainly contain species name and local temperature range.
7. In species four lines, second line contain 5 coefficients for high temperature range out of seven coefficients.(Highlighted in orange color)
8. In species four lines, third line contain two coefficients for high temperature range out of seven coefficients (Highlighted in orange color) And three coefficients for low temperature range out of seven coefficients (Highlighted in green color)
9. In species four lines, fourth line contain four coefficients for low temperature range out of seven coefficients (Highlighted in green color)

Here main objective is to extract this coefficient for particular species. After the extraction of species coefficient data calculation of thermodynamic property is done by using following formulas.

`C_p/R = a_1 + a_2 T + a_3 T^2 + a_4 T^3 + a_5 T^4`

`H/{RT} = a_1 + a_2 T/2 + a_3 T^2 /3 + a_4 T^3 /4 + a_5 T^4 /5 + a_6/T`

`S/R = a_1 lnT + a_2 T + a_3 T^2 /2 + a_4 T^3 /3 + a_5 T^4 /4 + a_7`

By using this formula ther thermodynamic properties such as Enthalpy, Entropy and Specific heat are calculated.

MATLAB Codes :-

For Calculating Thermodynamic property and extracting the coefficient data following codes are use in MATLAB and the logic of coding is explain after the codes.

To generate plot and calculate property Four functions program and one main program is coded by me.

1. Function program for Extracting coefficients from file.
2. Function program for calculating specific heat, Enthalpy and Entropy.
3. Function program for Calculating molecular weight.
4. Function program for generate figure and save the figure and display molecular weight as output in command window.
5. Main program to call all those function and generate output.

1. Function program for Extracting coefficients from file : -

Following function is use to get coefficient from the four lines which related to particular species and contains higher and lower coefficients.

function [low_coeff,high_coeff]= Coefficients(line1,line2, line3, line4)
% Three lines related to coefficient of Species are taken as input and lower and higher temperature coefficient give as output 

% Get higher coefficient
a1_high = str2num(line2(1:15)); 
a2_high = str2num(line2(16:30));
a3_high = str2num(line2(31:45));
a4_high = str2num(line2(46:60));
a5_high = str2num(line2(61:75));             
a6_high = str2num(line3(1:15));
a7_high = str2num(line3(16:30));

% Get lower coefficients
a1_low = str2num(line3(31:45));
a2_low = str2num(line3(46:60));   
a3_low = str2num(line3(61:75));
a4_low = str2num(line4(1:15)); 
a5_low = str2num(line4(16:30));
a6_low = str2num(line4(31:45));
a7_low = str2num(line4(46:60));

% Generate Vector of lower and Higher coefficients
low_coeff = [a1_low,a2_low,a3_low,a4_low,a5_low,a6_low,a7_low];
high_coeff = [a1_high,a2_high,a3_high,a4_high,a5_high,a6_high,a7_high];

end

% Save this file with name 'Coefficients.m'  important

2. Function program for calculating specific heat, Enthalpy and Entropy :-

This function program is created by using the equations to which specific heat, Enthalpy and Entropy can be measured by using coefficients. This equations are given in upper section of report.

function [Cp,H,S] = thermo_data_calculate(Temp,coefficient_value)
% Take temperature vector and coefficients vector as input and give calculated specific heat,Enthalpy and entropy

% Value of universal gas constant R = 8.314 J/(mol·K)
R= 8.314;

% Distribute the Input coefficient value vector into individual coefficients
a1 = coefficient_value(1);
a2 = coefficient_value(2);
a3 = coefficient_value(3);
a4 = coefficient_value(4);
a5 = coefficient_value(5);
a6 = coefficient_value(6);
a7 = coefficient_value(7);

% For every value of temperature vector calculate the Cp, H and S by using for loop.

for i=1:length(Temp)
T=Temp(i);
% Calculate Cp
Cp(i) = ( a1 + a2*T + a3*T^2 + a4*T^3 + a5*T^4) * R;
% Calculate Enthalpy
H(i) = ( a1 + a2*T/2 + a3*(T^2)/3 + a4*(T^3)/4 + a5*(T^4)/5 + a6/T) *R*T;
% Calculate Entropy
S(i) = ( a1*log(T) + a2*T + a3*(T^2)/2 + a4*(T^3)/3 + a5*(T^4)/4 + a7) *R;
end
end 

% Save this file with name 'thermo_data_calculate.m'

3. Function program for Calculating molecular weight :-

In this function values of atomic weights of atom is provided and by using this automic weight the program calculate molecular weight of input species.

function weight = molecular_weight(formula)
% Take Molecular formula as an input and give molecular weight as output. 

% Assign Automic weight for the elements which are use in data file.
H = 1.00784;
O = 15.999;
N = 14.0067;
C = 12.0107;
S = 32.065;
AR = 39.948;

% Calculate molecular weight
weight = 0;
for i=1:length(formula)
    if strcmp(formula(i),'H') % Compare 'i' element of formula with 'H'
        weight = weight + H ;
        W1 = H;
    elseif strcmp(formula(i),'O') % Compare 'i' element of formula with 'O'
        weight = weight + O ;
        W1 = O;
    elseif strcmp(formula(i),'N') % Compare 'i' element of formula with 'N'
        weight = weight + N ;
        W1 = N;
    elseif strcmp(formula(i),'C') % Compare 'i' element of formula with 'C'
        weight = weight + C ;
        W1 = C;
    elseif strcmp(formula,'CH2(S)')% Compare formula with 'CH2(S)'
        weight = C + 2*H + S;
        break; % We break this formula because there are only one species with this formula
    elseif strcmp(formula(i),'A') % Compare 'i' element formula with 'A' 
        weight = AR ;
        break; % We break this formula because there are only one species with this formula
    else 
        A = str2num(formula(i));% If there are number quantity in formula
        
        weight = weight + (A-1)*W1 ;
    end   

end

end

% Save this file with name 'molecular_weight.m'

4. Function program for generate figure and save the figure and display molecular weight as output in command window:-

This program generate output figure for input species and also save the figure in current directory and also print the molecular weight of species by calling molecular weight function.

function fig = figure_call_all(Input,i,lower_temp_boundary,upper_temp_boundary, Cp_low,Cp_high,H_low,H_high, S_low, S_high)    
% Take the data which are required to generate figures.

% Calculate molecular weight by calling molecular weight formula
Molecular_weight = molecular_weight(Input);

% Give command to generate figure in full screen
figure('units','normalized','outerposition',[0 0 1 1])

fig = figure(i);

% Generate subplot of 2 by 2
subplot(2,2,1) % First part of subplot is contain temerature v/s specific heat plot
plot(lower_temp_boundary,Cp_low,'linewidth',2,'color','k')
hold on
plot(upper_temp_boundary,Cp_high,'linewidth',2,'color','k')
xlabel('Temperature [K]')
ylabel('Specific Heat Cp  [ J/(mol-K) ]')
title(['Temperature V/S Specific Heat for ', Input])
   
subplot(2,2,2) % Second part of subplot is contain temerature v/s Enthalpy plot
plot(lower_temp_boundary,H_low,'linewidth',2,'color','r')
hold on
plot(upper_temp_boundary,H_high,'linewidth',2,'color','r')
xlabel('Temperature [K]')
ylabel('Enthalpy H  [ J/mol ]')
title(['Temperature V/S Enthalpy for ', Input])
      
subplot(2,2,3) % Third part of subplot is contain temerature v/s Entropy plot
plot(lower_temp_boundary,S_low,'linewidth',2,'color','b')
hold on
plot(upper_temp_boundary,S_high,'linewidth',2,'color','b')
xlabel('Temperature [K]')
ylabel('Entropy [ J/(mol-K) ]')
title(['Temperature V/S Entropy for ', Input])
   
M_W = num2str(Molecular_weight); % Convert molecular weight from number to string
str=['Molecular Mass of ', Input,' is ', M_W];
ax = subplot(2,2,4); % Fourth part of subplot is contain value of particular species molecular weight
text(0,0.5,str)
set(ax,'visible','off')
print(fig,[Input] ,'-djpeg') % Save figure in directory

fprintf('Molecular weight of %s is %f \n',Input, Molecular_weight) % Print molecular weight in command window
end

% Save this file with name 'figure_call_all.m'

5. Main program to call all those function and generate output :-

This is the main program by which all species figure and molecular weight is generated And we get the solution of our problem.

% Project : Parsing NASA thermodynamic data
clear all   % Clear the Workspace
close all   % Clear all figures and plot
clc         % Clear the Command Window

% Open file in MATLAB by using fopen command.
% 'r' is use for reading the file.(Operation)
Nasa_data = fopen('THERMO.dat','r');

fgetl(Nasa_data);   % Skip first header line 

% Read Second line Which contain Global Temperature Range                       
G_Temp_str = fgetl(Nasa_data);

% Split this string line into values of temperature by using strsplit command
G_Temp = strsplit(G_Temp_str,' ');  % Split by space in line

% This command split the line and save it in cell array as a string
% str2num command use to convert this string to numerical value.
% Values are taken from second index because first index value is a blank matrix
low_temp_global = str2num(G_Temp{2});
mid_temp_global = str2num(G_Temp{3});
high_temp_global = str2num(G_Temp{4});

for i=1:3
    % Skip next three lines which contain comments about THERMO.dat file
    % Which are not required
    fgetl(Nasa_data);
end

n=53; % Total number of species

for i=1:n
    line1=fgetl(Nasa_data); % Get first line of species
         
    line1_split = strsplit(line1,' '); % Split the line where space is located
    
    Input = line1_split{1}; % Take the name of species
    
    local_low_temp = str2num(line1(49:49+7)); % Calculate local low temperature
    local_hign_temp = str2num(line1(49+8:49+8+8)); % Calculate local high temperature
    local_mid_temp = str2num(line1(49+17:49+17+8)); % Calculate local middle temperature
        
    line2 = fgetl(Nasa_data); % Get second line of species
    line3 = fgetl(Nasa_data); % Get third line of species
    line4 = fgetl(Nasa_data); % Get four line of species
    
    % Generate coefficients by passing line1, line2, line3 and line4
    [low_coeff,high_coeff]= Coefficients(line1,line2, line3, line4);
    
    % Set lower and upper temperature boundaries
    lower_temp_boundary = local_low_temp : local_mid_temp;
    upper_temp_boundary = local_mid_temp : local_hign_temp;
    
    % Calculate Cp, H and S for lower and upper boundaries
    [Cp_high,H_high,S_high] = thermo_data_calculate(upper_temp_boundary,high_coeff);
    [Cp_low,H_low,S_low] = thermo_data_calculate(lower_temp_boundary,low_coeff);
    
    % Generate figure and molecular weight
    fig(i) = figure_call_all(Input,i,lower_temp_boundary,upper_temp_boundary, Cp_low,Cp_high,H_low,H_high, S_low, S_high);   
      
    
end

By run this final program we get desire output from it.

Program Explanation : -

• Main objective of this program is to extract coefficient data from a file which name is 'THERMO.dat' and find out Cp , H and S value from this data and generate there plot with respect to time.
• To do fopen and fgetl command is used for file parsing.
• First 'fopen' command is used to open the THERMO.dat file in matlab and give reading permission by using 'r' in command.
• After the opening of file 'fgetl' command is use for read lines of file. This command read file line by line with it's execution. For example when we execute it first time it read first line of file and after second time execution it read second line and simillary for n time execution it read n number line.
• In program first line of data file is not required for that we run the 'fgetl' command one time without assign any variable to it. After that second line contain global temperature range for that we run the 'fgetl' command second time and assign line value as a string to a varible.
• By using strsplit command we split the string with the location of space. and by using indexing and str2num command we assign value of second line to meaningful variables.
• Next three lines contain comments for that we skip this lines by calling 'fgetl' command and not assigning any variable to it.
• There are 53 species data contain by THERMO.dat file for that to extract all species data we use for loop to extract all species data.
• In every loop 'fgetl' command is used four times so that it extract the four lines which contain a species data. line1,line2,line3 and line4 are name given to those four lines.
• As known from problem statement section that line1 contain species name and local temperature range. By using approprite use of splitstr and indexing and str2num commands we get the species name and local temperature range.
• line2, line3 and line4 contain coefficient data. To get coefficient data a function is created. This function take this lines as input and give the two vectors as low temperature range and high temperature range vectors. Extraction is done by using indexing and then convert str format to numerical format by using str2num command.
• After the extraction of coefficients, temperature and coefficients are pass through second function which calculate the value of Cp, H and S for the temperature range by using above mentioned formula.
• After the calculation of thermodynamic property we pass this property to a new fuction which generate figure and save the figure in the directory with species name and also calculate the molecular weight by using a function created for calculating molecular weight.
• After generation of figure and calculation of molecular weight for one species then loop is run for second species and so on for other species.
• And finally figures and molecular weight is achived as an output.

Results :-

1. Screenshot of output folder

2. O2 Plot

3. N2 Plot

4. CO2 Plot

5. Molecular weight in Command window

6. Link of all result :

This link conatin all result figures.

https://drive.google.com/drive/folders/1NUkMlXlxY-CSrMx3gHDLz0V0GWqIQzA1?usp=sharing

Errors :

During programming syntax error and function passing errors are generated, those errors are removed by more study about commands and there syntax.

References :-

1. http://combustion.berkeley.edu/gri-mech/data/nasa_plnm.html
2. http://www.fkm.utm.my/~mazlan/?download=Comb%20-%20nasa.polynomial.pdf
3. http://combustion.berkeley.edu/gri-mech/version30/files30/thermo30.dat
4. MATLAB Documentation