Challenge

Project 1 - Parsing NASA thermodynamic data

AAkshat Nama11 Jun 2022

Aim:-

To write a code in MATLAB to parse the NASA thermodynamic data file and then calculate the thermodynamic properties of various gas species.

Objective:-

1. 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.

2. To calculate the molecular weight of each species and display it in the command window.

3. To plot the Cp, Enthalpy, and Entropy for the local temperature range (low temperature: high temperature) specific for each species.

4. To save the plots as images with appropriate names and dump them in separate folders for each species with a suitable name. All these should be generated automatically.

Theory:-

What is NASA Thermodynamic File:-

The NASA thermodynamics file, THERMO.dat, contains temperature ranges and 53 species, each of which contains 14 coefficients and temperature ranges, which we must parse to get Specific heat, Enthalpy, and Entropy using formula.

File Parsing:-

Syntactic analysis, also known as parsing syntax analysis, is the process of analysing a string of symbols that conforms to the norms of formal grammar, whether in natural language, computer languages, or data structures. Parsing is derived from the Latin word pars (orationis), which means part (of speech). The formal analysis by a computer of a phrase or other string of words into its constituents, resulting in a parse tree indicating their syntactic relationship to each other, which may also contain semantic and other information, is referred to as computational linguistics. For a syntactically ambiguous input, some parsing algorithms may generate a parse forest or a list of parse trees. Reading a data stream and constructing an in-memory model of the semantic content of that data is what parsing a file entails. This is intended to make it easier to do data transformations.

Specific Heat(Cp):-

The amount of heat per unit mass necessary to raise the temperature by one degree Celsius is known as specific heat. The relationship between heat and temperature change is represented below.

Enthalpy(H):-

A thermodynamic system's enthalpy is equal to the sum of the system's internal energy and the product of its pressure and volume. For operations at constant pressure, the heat absorbed or released equals the change in enthalpy in a system enclosed to avoid mass transfer.

Entropy(S):-

Entropy is the amount of thermal energy in a system per unit temperature that isn't being used to conduct beneficial work. Because work is generated by ordered molecular motion, entropy is also a measure of a system's molecular disorder or randomness.

Governing 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*lnT+a2*T+(a3*T^2)/2+(a4*T^3)/3+(a5*T^4)/4+a7`

Here R is the universal gas constant, T is the temperature, Cp is the specific heat(KJ), H is the enthalpy, and S is the entropy.

Body of the content(Code):-

% Parsing NASA Thermodynamic Data
clc 
clear all
close all

% Reading NASA Thermodynamic Data file
f1 = fopen('THERMO.dat','r')% Using 'fopen' command to open the file and using 'r' for reading the file 
get = fgetl(f1); % To read a line from file

% Reading Header Information
global_temp = fgetl(f1);% Obtaining the global temperatures

% Spliting 
temp = strsplit(global_temp);

% Using the str2num command to convert string cells to numbers
global_low_temp = str2num(temp{2});   % 2nd value of the temp is the global lowest temperature
global_medium_temp = str2num(temp{3});% 3rd value of the temp is the global medium temperature
global_high_temp = str2num(temp{4});  % 4th value of the temp is the global highest temperature

% Temperature Range
T = linspace(global_low_temp,global_high_temp,1000);

% Universal Gas Constant
R = 8.314;

% Skipping 3 comment lines 
for i = 1:3
    get = fgetl(f1);
end

% Reading species
for j = 1:53
    get = fgetl(f1);
    A = strsplit(get,' ');% Spliting into cell array
    species = (A{1});% To read the name of the species
    
    % Converting string to number 
    low_temp = str2num(A{end-3});
    midium_temp = str2num(A{end-2});
    high_temp = str2num(A{end-1});

    line1 = fgetl(f1);% Getting next line
    a = strfind(line1,'E');% Finding the string E
    
    % Finding First 7 Higher coefficients
    a1 = line1(1:a(1)+3);
    a1 = str2num(a1);
    a2 = line1(a(1)+4:a(2)+3);
    a2 = str2num(a2);
    a3 = line1(a(2)+4:a(3)+3);
    a3 = str2num(a3);
    a4 = line1(a(3)+4:a(4)+3);
    a4 = str2num(a4);
    a5 = line1(a(4)+4:a(5)+3);
    a5 = str2num(a5);
    
    line2 = fgetl(f1);
    b = strfind(line2,'E');
    
    a6 = line2(1:b(1)+3);
    a6 = str2num(a6);
    a7 = line2(b(1)+4:b(2)+3);
    a7 = str2num(a7);
    
    % Finding second 7 lower coefficients
    a8 = line2(b(2)+4:b(3)+3);
    a8 = str2num(a8);
    a9 = line2(b(3)+4:b(4)+3);
    a9 = str2num(a9);
    a10 = line2(b(4)+4:b(5)+3);
    a10 = str2num(a10);
    
    line3 = fgetl(f1);
    c = strfind(line3,'E');
    
    a11 = line3(1:c(1)+3);
    a11 = str2num(a11);
    a12 = line3(c(1)+4:c(2)+3);
    a12 = str2num(a12);
    a13 = line3(c(2)+4:c(3)+3);
    a13 = str2num(a13);
    a14 = line3(c(3)+4:c(4)+3);
    a14 = str2num(a14);
    
    % Calculate specific heat
    Cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,global_medium_temp);
    
    % Calculate Enthalpy
    H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,R,global_medium_temp);
    
    % Calculate Entropy
    S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,global_medium_temp);
    
    % Calculate molecular weight of the species
    W = molecular_weight(species);
    
    % Writing molecular weight file
    MW_file = fopen('Molecular weight.txt','w');
    fprintf(MW_file,'Molecular weight of %s is %d ',species,W);
    fclose(MW_file); % Closing the file
    
    % Saving in a folder
    cd('D:\SKILL-LYCN\Matlab\Practice\NASA Thermodynamic\New folder')
    mkdir(species);% Creating directory
    cd(species);
    
    % Plot1 Cp vs T
    figure(1)
    plot(T,Cp,'linewidth',2,'color','b');
    xlabel('Temperature(K)');
    ylabel('Specific Heat(KJ)');
    title(sprintf('Specific heat vs Temperature range of %s',species));
    grid on 
    saveas(figure(1),'Specific Heat.jpg');
    
    % Plot2 H vs T
    figure(2)
    plot(T,H,'linewidth',2,'color','g');
    xlabel('Temperature(K)');
    ylabel('Enthalpy(KJ/mol)');
    title(sprintf('Enthalpy vs Temperature range of %s',species));
    grid on 
    saveas(figure(2),'Enthalpy.jpg');
    
    % Plot3 S vs T
    figure(3)
    plot(T,S,'linewidth',2,'color','r');
    xlabel('Temperature(K)');
    ylabel('Entropy(KJ/mol-K)');
    title(sprintf('Entropy vs Temperature range of %s',species));
    grid on 
    saveas(figure(3),'Entropy.jpg');
    
    % Changing Directory
    cd('D:\SKILL-LYCN\Matlab\Practice\NASA Thermodynamic\New folder')
    
    % Molecular weight
    W = molecular_weight(species)
end

Explanation of code:-

1. Before starting the main code, I have used the commands. "clc" for clearing the command window, "clear all" for clearing all the previous variables, and "close all" for closing all the previous plots.

2. Then, using the reading commands 'fopen' and 'r,' read the text file Thermo.dat and open it.

3. Using 'fgetl' to read a line from a text file.

4. Split the temperature into a specific string cell using'strsplit' to get the temperature range from the Thermo.dat file.

5. 'str2num' was used to convert string cells to numbers.

6. Divided the temperature into a string and then transformed it to a number so that the NASA thermodynamics

could be solved independently using the temperature data.

7. The temperature range was given using 'linspace', which ranged from low to high. After that, add R, the universal gas constant.

8. Using "for & end loops" to skip the text file's following three comment lines.

9. Find the molecular weights for all 53 species by plotting specific heat(cp), enthalpy(H), and entropy(S) with temperature values.

10. Used to read 53 species in a loop.

11. I used 'strfind' to find a string.

12. Read line by line from the species to extract higher and lower temperature coefficients.

13. Using function codes for specific heat, enthalpy, and entropy.

14. Using 'mkdir' with a location URL to create a directory.

15. Using the plot command, plot all three plots for each species. Then there's the Molecular weight calculation.

Functions Used:-

Function for calculating Specific Heat:-

% Function for calculating Specific Heat

function Cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,global_medium_temp)

if (T <= global_medium_temp) 
    Cp = R*(a8 + (a9*T) + (a10*(T).^2) + (a11*(T).^3) + (a12*(T).^4));% With the help of given formula
else
    Cp = R*(a1 + (a2*T) + (a3*(T).^2) + (a4*(T).^3) + (a5*(T).^4));
end
end

Explanation:-

Will calculate the Specific heat by generating a Specific heat function and assigning a formula coefficient from 14 coefficients, including the Gas constant, Temperature, and Mid-temperature.To specify a lower temperature coefficient in the calculation, we'll use a 'if' loop to construct a condition with a temperature less than and equal to mid-temperature. Then, in the formula, we'll define the coefficient of the Higher temperature coefficient using the 'else' condition. Then end the loop and use it to Calculate the Specific Heat function.

Function for calculating Enthalpy:-

% Function for calculating Enthalpy

function H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,R,global_medium_temp)
if (T <= global_medium_temp)
    H = R*T.*(a8 + ((a9*T)/2) + ((a10*(T).^2)/3) + ((a11*(T).^3)/4) + ((a12*(T).^4)/5) + (a13./T));% With the help of given formula
else
    H = R*T.*(a1 + ((a2*T)/2) + ((a3*(T).^2)/3) + ((a4*(T).^3)/4) + ((a5*(T).^4)/5) + (a6./T));
end
end

Explanation:-

Will calculate the Enthalpy by generating an Enthalpy function and assigning a formula coefficient from 14 coefficients, including the Gas constant, Temperature, and Mid-temperature.To specify a lower temperature coefficient in the calculation, we'll use a 'if' loop to construct a condition with a temperature less than and equal to mid-temperature. Then, in the formula, we'll define the coefficient of the Higher temperature coefficient using the 'else' condition. Then end the loop and use it to Calculate the Enthalpy function.

Function for Entropy:-

% Function for Entropy

function S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,global_medium_temp)

if (T <= global_medium_temp)
    S = R*(a8.*log(T) + (a9*T) + ((a10*(T).^2)/2) + ((a11*(T).^3)/3) + ((a12*(T).^4)/4) + a14);% With the help of given formula
else
    S = R*(a1.*log(T) + (a2*T) + ((a3*(T).^2)/2) + ((a4*(T).^3)/3) + ((a5*(T).^4)/4) + a7);
end
end

Explanation:-

Will calculate the Entropy by generating an Entropy function and assigning a formula coefficient from 14 coefficients, including the Gas constant, Temperature, and Mid-temperature.To specify a lower temperature coefficient in the calculation, we'll use a 'if' loop to construct a condition with a temperature less than and equal to mid-temperature. Then, in the formula, we'll define the coefficient of the Higher temperature coefficient using the 'else' condition. Then end the loop and use it to Calculate the Entropy function.

Function for calculating Molecular weight:-

% Function for calculating Molecular weight

function W = molecular_weight(species)

% Defining atomic name of the species
% Hydrogen, Carbon, Oxygen, Nitrogen, Argon
atomic_name = ['H','C','O','N','Ar'];

% Defining atomic mass of the main species
atomic_weight = [1.008 12.011 15.999 14.007 39.948];
 
W = 0;% Starting molecular weight
for i = 1:length(species)% Loop runs to the length of the species
    for j = 1:length(atomic_name)% Loop runs to the length of the atomic names
        if strcmp(species(i), atomic_name(j))% By using 'strcmp' function we can compare two strings
            W = W + atomic_weight(j);
            position = j;
        end
    end
    n = str2num(species);% Now we'll see if there are more elements in the species, and if so, we'll increase the species' atomic weight 
    
    if n>1
        W = W + atomic_weight(position)*(n-1);
    end 
end

fprintf('Molecular Weight of %s : ',species);
fprintf('%f',W);
disp(' ')
end

Explanation:-

1. We begin by creating a Molecular weight function. W, and putting species input into it

2. Then, in an array atomic_name, define the major element of all 53 species, which are 'H', 'C', 'O', 'N', and 'Ar'.

3. Then, using the primary elements listed above, we'll create an atomic weight array with the values 1.008, 12.011, 15.999, 14.007, and 39.948.

4. Then we'll set the molecular beginning range W to zero and establish a for loop condition that runs the length of the species, followed by a nested for loop that runs the length of the string atomic name.

5. Then, using the 'if' condition, we'll compare the atomic_name strings and determine molecular masses using the 'strcmp' function. If both strings are found to be the same, the molecular weight will be assigned.

6. Now we'll use the 'str2num' command to convert the string to a number and assign it to n. Then we use 'if' loop to determine whether the species has more than one element; if there are more than one element, we will increase the atomic weight by multiplying the position by n-1 and adding it to the atomic_weight array.

7. The Molecular weight will then be printed using the 'fprintf' command and displayed using the 'disp' command.

Errors faced:-

1. Previously I am using curve brackets instead of curly brackets.