Project 1 - Parsing NASA thermodynamic data

AIM:

To determine the thermodynamic properties of various species using NASA thermodynamic data.

GOVERNING EQUATIONS:

`c_p/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 = a1lnT + a2*T + a3*T^2/2 + a4*T^3/3 + a5*T^4/4 + a7`

where,

Cp is specific heat in [kJ/kg-K]

H is enthalpy in [kJ/mol]

S is entropy in [kJ/mol-K]

T is temperature in Kelvin[K]

R is universal gas constant

R = 8.314 [kJ/mol-K]

OBJECTIVES:

1. To extract the high temperature and low temperature coefficients from the NASA thermodynamic data file.
2. To determine the molecular weight of each species.
3. To determine the enthalpy, entropy and specific heat of each species.
4. To plot the curve of the above properties against temperature.
5. To save the plots as images in a separate folder named as the species name.

MATLAB CODE:

1. Main script

clear all
close all
clc

% Universal Gas Constant( kJ/mol-K)
R = 8.314; 

% Opening NASA file
file = fopen('THERMO.dat','r');

% Reading Header Information
fgetl(file);
global_value = fgetl(file);
A = strsplit(global_value, '  ');

% Extracting Global Temperature Values
global_low_temp = str2double(A{2});
global_mid_temp = str2double(A{3});
global_high_temp = str2double(A{4});

% Skipping next 3 lines
fgetl(file);
fgetl(file);
fgetl(file);

for i = 1:53
    % Extracting species name
    
    line1 = fgetl(file);
    b = strsplit(line1,'  ');
    species_name = b{1};
    n = length(b);
    
    % Calculating Molecular Weight
    
    Mass = MolecularWeight(species_name);
    disp(['Molecular Weight of ', species_name,' is ',num2str(Mass)]);
    
    % Extracting Local temperature
    
    local_low_temp = str2double(b{n - 3});
    local_high_temp = str2double(b{n - 2});
    local_mid_temp = str2double(b{n - 1});
    
    %Extracting Coefficients
    line2 = fgetl(file);
    a = strfind(line2,'E');
    
    % High Temperature Coefficients
    
    a1 = str2double(line2(1:a(1) + 3));
    a2 = str2double(line2(a(1) + 4 : a(2) +3));
    a3 = str2double(line2(a(2) + 4 : a(3) +3));
    a4 = str2double(line2(a(3) + 4 : a(4) +3));
    a5 = str2double(line2(a(4) + 4 : a(5) +3));
    
    line3 = fgetl(file);
    a = strfind(line3,'E');
    
    a6 = str2double(line3(1:a(1) + 3));
    a7 = str2double(line3(a(1) + 4 : a(2) +3));
    
    % Low Temperature Coefficients
    
    a8 = str2double(line3(a(2) + 4 : a(3) +3));
    a9 = str2double(line3(a(3) + 4 : a(4) +3));
    a10 = str2double(line3(a(4) + 4 : a(5) +3));
    
    line4 = fgetl(file);
    a = strfind(line4,'E');
    a11 = str2double(line4(1:a(1) + 3));
    a12 = str2double(line4(a(1) + 4 : a(2) +3));
    a13 = str2double(line4(a(2) + 4 : a(3) +3));
    a14 = str2double(line4(a(3) + 4 : a(4) +3));
    
    % Creating Temperature Ranges
    T_low = linspace(local_low_temp, local_mid_temp, 200);
    T_high = linspace(local_mid_temp, local_high_temp, 300);
    
    % For High Temperature Ranges
    Cp_high = (a1 + a2 * T_high + a3 * T_high.^2 + a4 * T_high.^3 + a5 * T_high.^4) * R; 
    H_high = (a1 + a2 * T_high / 2 + a3 * T_high.^2 / 3 + a4 * T_high.^3 / 4 + a5 * T_high.^4 / 5 + a6 ./ T_high) * R .* T_high;
    S_high = (a1 * log(T_high) + a2 * T_high + a3 * T_high.^2 / 2 + a4 * T_high.^3 / 3 + a5 * T_high.^4 / 4 + a7) * R;
    
    % For Low Temperature Ranges
    Cp_low = (a8 + a9 * T_low + a10 * T_low.^2 + a11 * T_low.^3 + a12 * T_low.^4) * R; 
    H_low = (a8 + a9 * T_low / 2 + a10 * T_low.^2 / 3 + a11 * T_low.^3 / 4 + a12 * T_low.^4 / 5 + a13 ./ T_low) * R .* T_low;
    S_low = (a8 * log(T_low) + a9 * T_low + a10 * T_low.^2 / 2 + a11 * T_low.^3 / 3 + a12 * T_low.^4 / 4 + a14) * R;
    
    % Plotting
    % Specific Heat
    figure(1)
    plot(T_low, Cp_low, T_high, Cp_high, 'color', 'r', 'linewidth', 3);
    title(['Temperature vs Specific Heat of ', species_name], 'fontsize', 16);
    xlabel('Temperature [K]', 'fontsize', 14);
    ylabel('Heat Capacity [kJ/kg-K]', 'fontsize', 14);
    grid on;
    
    % Enthalpy
    figure(2)
    plot(T_low, H_low, T_high, H_high, 'color', 'b', 'linewidth', 3);
    title(['Temperature vs Enthalpy of ', species_name], 'fontsize', 16);
    xlabel('Temperature [K]', 'fontsize', 14);
    ylabel('Enthalpy [kJ/mol]', 'fontsize', 14);
    grid on;
    
    % Entropy
    figure(3)
    plot(T_low, S_low, T_high, S_high, 'color', 'k', 'linewidth', 3);
    title(['Temperature vs Entropy of ', species_name], 'fontsize', 16);
    xlabel('Temperature [K]', 'fontsize', 14);
    ylabel('Entropy [kJ/mol-K]', 'fontsize', 14);
    grid on;
    
    % Saving Files
    original_path = pwd;
    mkdir('D:\Skill-Lync\MATLAB\File Parsing',species_name)    % Creating directory
    cd (['D:\Skill-Lync\MATLAB\File Parsing\',species_name])
    
    saveas(1,'Specific Heat.jpeg');
    saveas(2,'Enthalpy.jpeg');
    saveas(3,'Entropy.jpeg');
    
    mole = fopen('Weight.txt', 'w');
    fprintf(mole, 'Molecular Weight of  %s is %d', species_name, Mass);
    fclose(mole);
    
    cd (original_path)
    
end

% Closing NASA file
fclose(file);

fopen command is used to load the NASA thermodynamic data file.

fgetl is used to read a particular line.

strsplit command is used to split the line according to some given arguments.

strfind is used to find the position of a particular character.

str2double is used to convert the string into numeric data type.

mkdir commands create a new directory.

cd command changes the directory.

pwd commands return the path of current directory.

saveas command is used to save the plots as an image file in the current folder.

fclose closes the specified file.

I have created to temperature ranges, one for high ranges and one for low ranges.

2. Molecular Weight File

function fval = MolecularWeight(species_name)
 % Calculating Molecular Weight of the species
    mass_data = struct('C', 12, 'O', 16, 'H', 1, 'N', 14, 'A', 40);
    Mass = 0;
    j = 1;
    while j <= length(species_name)
        if species_name(j)=='C' || species_name(j)=='O' || species_name(j)=='H' || species_name(j)=='N' || species_name(j)=='A'
            if j < length(species_name)
                num = str2double(species_name(j + 1));
            else
                num = 1;
            end
            
            if num>1
                Mass = Mass + mass_data.(species_name(j)) * num;
                j = j+1;
            else
                Mass = Mass + mass_data.(species_name(j));
            end
        end
        j = j + 1;
    end
    fval = Mass;
end

I have used struct command to create a structure containing the first letter of elements as its name and molar mass as its value.

Before entering into the inner if statement I have checked whether the particular character is the first letter of the elements or not so that it will not result in an error while searching the value in the structure.

I have checked whether the next character is a number in the species name of an element so that the program can calculate the exact mass of the species without showing an error.

This program can only calculate the mass of the molecule containing the five elements mentioned in the file.

OUTPUT:

1. Screenshot of the window where all the folders are saved

2. Plots for O2

a) Specific Heat

b) Enthalpy

c) Entropy

3. Plots for N2

a) Specific Heat

b) Enthalpy

c) Entropy

4. Plots for Co2

a) Specific Heat

b) Enthalpy

c) Entropy