Challenge

Project 1 - Parsing NASA thermodynamic data

SSuraj Dhavaleshwar11 Nov 2021

NASA Thermodynamic data file:- https://drive.google.com/file/d/1twwGJA-XQfAnk2DXcKVb-wduU2i7yd4t/view?usp=sharing

PARSING NASA THERMODYNAMIC DATA

AIM: To write a function to extract coefficients from a data file and calculate the specific heats, enthalpy, entropy and molecular weight for all the species in the given data file

OBJECTIVE:

To read THERMO.dat file
To code a function which extracts the coefficients and calculates the enthalpy, entropy, and specific heats
To calculate the molecular weight and display in the command window
Plot the Cp, Enthalpy and Entropy for the local temperature range specific for each species
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

FORMULAE:

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

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

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

Where,

Cp = Specific heat at constant pressure, J/(kg.K)

H = Enthalpy, J

S = Entropy, J/K

T = Temperature, degrees

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

Specific heat(Cp): It is the amount of heat that must be added to one unit of mass of the substance in order to cause an increase of one unit in temperature. The SI unit of specific heat is joule per kelvin per kilogram, J/(kg.K)

Enthalpy(H): It is defined as the sum of the system's internal energy and the product of its pressure and volume. The SI unit of enthalpy is joule, J

Entropy(S): It is the measure of a system’s thermal energy per unit temperature that is unavailable for doing useful work. Because work is obtained from ordered molecular motion, the amount of entropy is also a measure of the molecular disorder, or randomness, of a system. The SI unit of entropy is joule per kelvin, J/K

BODY OF CONTENT:

Programming Language: MATLAB

Code:

NASA Thermodynamic data file reading

% NASA Thermodynamic data
clear all
close all
clc

F1 = fopen('THERMO.dat','r');       % opening existing file in read mode
file_name = fgetl(F1);              % reading first line 
A = fscanf(F1,'%f');                % file reading           
A = A';                             % transpose of column vector
% Assigning the data to variables
global_low_temp = A(1);
global_mid_temp = A(2);
global_high_temp = A(3);
    
    
for i = 1:3
    comments = fgetl(F1);           % reading next three lines using fgetl command in for loop
end


for j = 1:53                 
    l1 = fgetl(F1);                 % assigning first line to variable l1 
    L = strsplit(l1,' ');           % splitting the line with space as delimiter
    l = length(L);
    species_name = L{1};
    local_low_temp = str2double(L(l-3));
    local_mid_temp = str2double(L(l-2));
    local_high_temp = str2double(L(l-1));
    
    l2 = fgetl(F1);                  % assigning second line to variable l2
    x = strfind(l2, 'E');            % finding positions of E 
    a1 = str2double(l2(1:x(1)+3));
    a2 = str2double(l2((x(1)+4):(x(2)+3)));
    a3 = str2double(l2((x(2)+4):(x(3)+3)));
    a4 = str2double(l2((x(3)+4):(x(4)+3)));
    a5 = str2double(l2((x(4)+4):(x(5)+3)));
    
    l3 = fgetl(F1);                   % assigning third line to variable l3
    y = strfind(l3,'E');              % finding positions of E
    a6 = str2double(l3(1:y(1)+3));
    a7 = str2double(l3((y(1)+4):(y(2)+3)));
    a8 = str2double(l3((y(2)+4):(y(3)+3)));
    a9 = str2double(l3((y(3)+4):(y(4)+3)));
   a10 = str2double(l3((y(4)+4):(y(5)+3)));
    
    l4 = fgetl(F1);                   % assigning fourth line to variable l4
    z = strfind(l4,'E');              % finding positions of E
    a11 = str2double(l4(1:z(1)+3));
    a12 = str2double(l4((z(1)+4):(z(2)+3)));
    a13 = str2double(l4((z(2)+4):(z(3)+3)));
    a14 = str2double(l4((z(3)+4):(z(4)+3)));
    
    
    % Calculation of Specific heat(Cp), Enthalpy(H), and Entropy(S)
    R = 8.314;                                              % universal gas constant,(J/mol.K)
    temp = linspace(local_low_temp,local_high_temp,1000);   % low temperature range

    Cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,R,global_mid_temp);
    H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,temp,R,global_mid_temp);
    S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,R,global_mid_temp);
    
    % calculation of Molecular weight
    
    M_W = molecular_weight(species_name);
    %M = fopen('molecular weight.txt','a+'); % creating file 'molecular weight' 
    fprintf('Molecular weight of %s = %f un', species_name, M_W)
    %fclose(M);
    
    
    
    % Plotting
    % Folder for saving plots
    maindir = pwd;      % displaying current working directory and assigning it to variable maindir
    mkdir(species_name);% make new directory
    cd(species_name);   % change current working directory to newly created one
    
    % specific heat(Cp) vs Temperature(temp)
    figure(1)
    plot(temp,Cp,'linewidth',2,'color','r');
    xlabel('Temperature[K]');
    ylabel('Specific heat[kJ]');
    title(sprintf('Specific heat vs Temperature of %s',species_name));
    saveas(figure(1),'Specific heat vs Temperature.jpg');
    
    % Enthalpy(H) vs Temperature(temp)
    figure(2)
    plot(temp,H,'linewidth',2,'color','k');
    xlabel('Temperature,[K]');
    ylabel('Enthalpy,[kJ/mol]');
    title(sprintf('Enthalpy vs Temperature of %s',species_name));
    saveas(figure(2),'Enthalpy vs Temperature.jpg');
    
    % Entropy(S) vs Temperature(temp)
    figure(3)
    plot(temp,S,'linewidth',2,'color','g');
    xlabel('Temperature,[K]');
    ylabel('Entropy,[kJ/(mol.k)]');
    title(sprintf('Entropy vs Temperature of %s',species_name));
    saveas(figure(3),'Entropy vs Temperature.jpg');
    
    cd(maindir);    % change directory to the old one
end

Function to calculate Specific heat(Cp)

% Function template to calculate Specific heat, Cp

% Specific heat, Cp

function Cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,R,global_mid_temp)
% Comparing temperature with medium global temperature to decide if low or
% high temperature coefficients should be used
if temp <= global_mid_temp
    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

Function to calculate Enthalpy(H)

% Function template to calculate Entalpy,H

% Enthalpy,H

function H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,temp,R,global_mid_temp)
% Comparing temperature with medium global temperature to decide if low or
% high temperature coefficients should be used
if temp <= global_mid_temp
    H = R*temp.*(a8+((a9*temp)/2)+((a10*(temp).^2)/3)+((a11*(temp).^3)/4)+((a12*(temp).^4)/5)+((a13./temp)));
else
    H = R*temp.*(a1+((a2*temp)/2)+((a3*(temp).^2)/3)+((a4*(temp).^3)/4)+((a5*(temp).^4)/5)+((a6./temp)));
end
end

Function to calculate Entropy(S)

% Function template to calculate Entropy,S

% Entropy,S

function S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,R,global_mid_temp)
% Comparing temperature with medium global temperature to decide if low or
% high temperature coefficients should be used
if temp <= global_mid_temp
    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

Function to calculate Molecular weight(u)

% Function to calculate Molecular weight

function M_W = molecular_weight(species_name)

element = {'H','C','N','O','AR'};  
atomic_weight = [1.00784, 12.0107, 14.0067, 15.999, 39.948];    % Atomic weights
M_W = 0;    % initializing molecular weight to 0

% Molecular weight calculation

for i = 1:length(species_name)
    for j = 1: length(atomic_weight)
        if strcmp(species_name(i), element(j))
            M_W = M_W +atomic_weight(j);
            u = j;
        end
    end
    n = str2double(species_name(i));
    if n>1
        M_W = M_W + (atomic_weight(u)*(n-1));
    end
end
end

Explaination:

Initially the workspace, plot windows and command screen are cleared using clear all, close all and clc respectively
The Thermo.dat file is being opened for reading using fopen in read mode and assigned to variable F1
First line is read using fgetl and assigned to variable file_name
Next three lines contains numbers, so these three lines are read using fscanf and stored in variable A and then the row vector A is transposed to get column vector
This column vector has three rows each having a number, these numbers are global temperature range and so are assigned to variables global_low_temp, global_mid_temo, and global_high_temp respectivey
Next three lines are read using fgetl command in for loop and stored in variable comments
Another for loop is used to read the remaining data in the file, in this loop the data is read in segments where each segment consists of four lines
From first line species name and local temperature range is extracted and stored and from the remaining three lines fourteen numbers are extracted which are coefficients used to calculate the specific heat, enhalpy and entropy and are stored in variables a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13, and a14 respectively
This is iterated for 53 times as there are 53 segments containing information of 53 species
Once the coefficients are extracted, the functions for calculating specific heat, enthalpy, and entropy is written which takes the 14 coefficients, local temperature range, R(universal gas constant) and global mid temperature as inputs and calculates the specific heat, enhalpy and entropy
The function for calculating the molecular weight is also written which takes species name as input, in this function the elements and there atomic weight are stored in two arrays namely element and atomic weight
The output variable M_W which is going to store the output value of function is initialized to zero
Then the nested for loop is used where the string comparision is done and molecular weight is calculated and if the species name contain more than one element then the molecular weight is updated using the if condition and final molecular weight is stored in output variable
Once all functions are written they are called inside main code
The molecular weight of each species is displayed in the command window
Plots of specific heat, enthalpy and entropy is ploted for the local temperature range and the plots of each species are stored in a seperate folder using mkdir() and cd() commands which are used to make new directory and change working directory respectively

ERRORS:

This is a case sensitive error, as I used Temp instead of temp, so it is displaying Unrecognised function or variable, after replacing Temp with temp this error is solved

In the below figure we can see the command window output, where the molecular weight of each species is displayed, but the format is wrong and it is repeated. This was due to incorrect loop condition, after changing the loop condition this problem was solved

OUTPUT:

Command window output showing Molecular weight of each species from the NASA Thermodynamic data file

Molecular weight of O = 15.999000 u
Molecular weight of O2 = 31.998000 u
Molecular weight of H = 1.007840 u
Molecular weight of H2 = 2.015680 u
Molecular weight of OH = 17.006840 u
Molecular weight of H2O = 18.014680 u
Molecular weight of HO2 = 33.005840 u
Molecular weight of H2O2 = 34.013680 u
Molecular weight of C = 12.010700 u
Molecular weight of CH = 13.018540 u
Molecular weight of CH2 = 14.026380 u
Molecular weight of CH2(S) = 14.026380 u
Molecular weight of CH3 = 15.034220 u
Molecular weight of CH4 = 16.042060 u
Molecular weight of CO = 28.009700 u
Molecular weight of CO2 = 44.008700 u
Molecular weight of HCO = 29.017540 u
Molecular weight of CH2O = 30.025380 u
Molecular weight of CH2OH = 31.033220 u
Molecular weight of CH3O = 31.033220 u
Molecular weight of CH3OH = 32.041060 u
Molecular weight of C2H = 25.029240 u
Molecular weight of C2H2 = 26.037080 u
Molecular weight of C2H3 = 27.044920 u
Molecular weight of C2H4 = 28.052760 u
Molecular weight of C2H5 = 29.060600 u
Molecular weight of C2H6 = 30.068440 u
Molecular weight of CH2CO = 42.036080 u
Molecular weight of HCCO = 41.028240 u
Molecular weight of HCCOH = 42.036080 u
Molecular weight of H2CN = 28.033080 u
Molecular weight of HCN = 27.025240 u
Molecular weight of HNO = 31.013540 u
Molecular weight of N = 14.006700 u
Molecular weight of NNH = 29.021240 u
Molecular weight of N2O = 44.012400 u
Molecular weight of NH = 15.014540 u
Molecular weight of NH2 = 16.022380 u
Molecular weight of NH3 = 17.030220 u
Molecular weight of NO = 30.005700 u
Molecular weight of NO2 = 46.004700 u
Molecular weight of HCNO = 43.024240 u
Molecular weight of HOCN = 43.024240 u
Molecular weight of HNCO = 43.024240 u
Molecular weight of NCO = 42.016400 u
Molecular weight of CN = 26.017400 u
Molecular weight of HCNN = 41.031940 u
Molecular weight of N2 = 28.013400 u
Molecular weight of AR = 0.000000 u
Molecular weight of C3H8 = 44.094820 u
Molecular weight of C3H7 = 43.086980 u
Molecular weight of CH3CHO = 44.051760 u
Molecular weight of CH2CHO = 43.043920 u

O2 plots:

1. Specific heat vs Temperature

2. Enthalpy vs Temperature

3. Entropy vs Temperature

N2 plots:

1. Specific heat vs Temperature

2. Enthalpy vs Temperature

3. Entropy vs Temperature

CO2 plots:

1. Specific heat vs Temperature

2. Enthalpy vs Temperature

3. Entropy vs Temperature

Window screenshot where all folders are saved

CONCLUSION:

The functions to calculate the specific heat, enthalpy, and entropy has been written and the NASA thermodynamic data file is read and the required data is extracted from the data file and the specific heat, enthalpy and entropy are calculated. The function to calculate the molecular weight is also written and called inside the main code to calculate the molecular weight and displayed in the command window. The plots are plotted and are stored in seperate folders for each species

REFERENCES:

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.