Project 1 - Parsing NASA thermodynamic data

Aim :- To write matlab code to obtain thermodynamic properties of different chemical 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 NASA thermodynamic data file "Thermo.dat".

The NASA polynomials have the form as shown below:

Where R is the universal gas constant, T is the temperature

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 suitable name automatically.

Function for Enthalpy plot:-

function[h] = enthalpy_file(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,R,temp,local_temp_mid)
        if temp <= local_temp_mid
           h = R*temp.*(b1+ (b2.*(temp./2))+ (b3.*(((temp).^2)./3))+ (b4.*(((temp).^3)./4))+ (b5.*(((temp).^4)./5)) + b6./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 for Entropy plot:-

function[s] = entropy_file(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,R,temp,local_temp_mid)
        if temp<= local_temp_mid
           s = R.*(b1.*log(temp) + (b2.*temp) + (b3.*(((temp).^2)./2)) + (b4.*(((temp).^3)./3)) + (b5.*(((temp).^4)./4)) + b7);
        else
           s = R.*(a1.*log(temp) + (a2.*temp) + (a3.*(((temp).^2)./2)) + (a4.*(((temp).^3)./3)) + (a5.*(((temp).^4)./4)) + a7);
        end
end

Function for Specific Heat plot:-

function[c_p] = sp_heat_file(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,temp,local_temp_mid)
        if temp <= local_temp_mid
           c_p = (b1+ (b2.*temp)+ (b3.*((temp).^2))+ (b4.*((temp).^3))+ (b5.*((temp).^4)));
        else
           c_p = (a1+ (a2.*temp)+ (a3.*((temp).^2))+ (a4.*((temp).^3))+ (a5.*((temp).^4)));
        end
end

Function for Molecular Weight plot:-

function[Mw] = molecular_weight_file(species_name)

    elements = ['H','C','N','O','A']
    atomic_weight = [1,12,14,16,40]

    Mw = 0;
    for i = 1:length(species_name)
        for j = 1:length(elements)
            if strcmp(species_name(i), elements(j)) 
               Mw = Mw + atomic_weight(j);
               a = j;
            end
        end
        n = str2double(species_name(i));
        if n>1
           Mw = Mw+ atomic_weight(a)*(n-1);
        end
    end
end

Main Code:-

clear all;
close all;
clc;

% Open the file for reading
f1 = fopen('THERMO.DAT','r');
line1 = fgetl(f1);
disp(line1);

line2 = fgetl(f1);
disp(line2);

% Read the data from the file
%data = fscanf(f1,'%f');

%temp = fgetl(f1);

% splitting 
s = strsplit(line2,' ');
R = 8.314;

% converting string cells into numbers
temp_low = str2num(s{2});
temp_mid = str2num(s{3});
temp_high = str2num(s{4});

%temperature Range
% skipping comment lines

for i = 1:3
    data = fgetl(f1);
    disp(data)
end

for j = 1:53
    Line1= fgetl(f1);
    species_desc = strsplit(Line1,' ');
    species_name = species_desc{1};

    a =flip(species_desc);
    local_temp_low = str2num(s{4});
    local_temp_mid = str2num(s{3});
    local_temp_high = str2num(s{2});

    temp = linspace(local_temp_low, local_temp_high,500);

    Line2 = fgetl(f1);
    a = findstr(Line2,'E');
    Line3 = fgetl(f1);
    b = findstr(Line3,'E');
    Line4 = fgetl(f1);
    c = findstr(Line4,'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));
    a6 = str2double(Line3(1:b(1)+3));
    a7 = str2double(Line3(b(1)+4:b(2)+3));
    
    %Low temeprature coefficients
    b1 = str2double(Line3(b(2)+4:b(3)+3));
    b2 = str2double(Line3(b(3)+4:b(4)+3));
    b3 = str2double(Line3(b(4)+4:b(5)+3));
    b4 = str2double(Line4(1:c(1)+3));
    b5 = str2double(Line4(c(1)+4:c(2)+3));
    b6 = str2double(Line4(c(2)+4:c(3)+3));
    b7 = str2double(Line4(c(3)+4:c(4)+3));
    
    Enthalpy = enthalpy_file(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,R,temp,local_temp_mid);
    Entropy = entropy_file(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,R,temp,local_temp_mid);
    Molecular_weight = molecular_weight_file(species_name);
    Specific_heat = sp_heat_file(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,temp,local_temp_mid);

    figure(1)
    plot(temp,Enthalpy,'linewidth',2,'color','b')
    title('temp vs enthalpy',species_name)
    xlabel('temperature')
    ylabel('Enthalpy')
    grid on
        
    figure(2)
    plot(temp,Entropy,'linewidth',2, 'color','r')
    title('temp vs entropy',species_name)
    xlabel('temp')
    ylabel('Entropy')
    grid on
      
    
    figure(3)
    plot(temp,Specific_heat,'linewidth',2,'color','y')
    title('temp vs specific heat',species_name)
    xlabel('temp')
    ylabel('specific heat')
    grid on
    
    mkdir(species_name);
    cd(species_name);
     
    saveas(1,'enthalpy.jpg')
    saveas(2,'entropy.jpg')
    saveas(3,'specific_heat.jpg')
    cd  ..
    disp(['molecular weight of ',species_name]);
    disp(Molecular_weight)
end

% Close the file
fclose(f1);

Explanations :-

• First three lines in main code is used for clear the command window,workspace.
• Using 'fopen' and 'r' for the opening thermo.dat file and r standa for reading the file.
• To obtain temperature range from the file use strsplit command that are split temperature into the string cell.
• Adding universal gas constant R.
• Used 'str2num' for converting string cell to numbers.
• Used 'findstr' to find string.
• For each species, the code reads several lines of data that contain coefficients for calculating enthalpy, entropy, and specific heat. These coefficients are extracted by parsing the lines of data using the 'findstr' and 'str2double' functions.
• Once the coefficients are obtained, the code calculates the values of enthalpy, entropy, specific heat, and molecular weight for a range of temperatures using functions named 'enthalpy_file', 'entropy_file', 'sp_heat_file', and 'molecular_weight_file', respectively.
• Use of mkdir command for making a directory and saveas for saving images for thermo property.

Output:-

Directory -

1. Plots for CO2

a. Variation of Specific heat with the temperature of O2

b. Variation of Entropy with the temperature of CO2

c. Variation of Enthalpy with the temperature of CO2

2. Plots for N2

a. Variation of Specific heat with the temperature of N2

b. Variation of Entropy with the temperature of N2

c. Variation of Enthalpy with the temperature of N2

Link for the folders :- https://drive.google.com/drive/folders/1LFhBSj8dPmy5EJdBQHMkRHDmgpfuMnD6?usp=share_link

Conclusion :

Using MATLAB code, the necessary coefficients required to evaluate various thermodynamic parameters for different species are extracted from NASA thermodynamic data file. The specific heat, enthalpy, and entropy for all the species are calculated and plotted against the local temperature range. The molecular weight is calculated for all the species.