Parsing NASA thermodynamic data using MATLAB

Aim : Parsing NASA Thermodynamic Data in MATLAB

Data Parsing : It refers to the process of converting information represented in one form into another form to make it easier to work with. It is done by first partially analysing the data, understanding its logic & structure , then representing that structure in some code.

Governing Equations : Properties like Specific Heat(C_p), Enthalpy(H) & Entropy(S) of a particular element or compound changes with respect to change in temperature.

The following equations represent that below,

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

`H/(R*T) = 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`

R is the universal gas constant & T is the absolute thermodynamic temperature.

The NASA thermodynamic file contains these (a1-a7) coefficients for low temperature & high temperature range.

The data file link is https://drive.google.com/file/d/1AikCHDlz9s_fu8whE81qee_GD4AyYpCP/view

Objectives :

i. Write a function that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the data file.

ii. Calculate the molecular weight of each species and display it in the command window.

iii. Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species.

iv. Save the plots as images with appropriate names and dump them in separate folders for each species with suitable name(automatically).

Explanation of Code :

1. First, the 'THERMO.dat' file is loaded , then first line is skipped using fgetl, From second line value of Global Temperature are extracted & again next 3 lines are skipped.
2. Now several matrix are predefined for making code more efficient, To get the element names line1 gets the first line required contanning the name of compounds & speecific local temperature range.
3. Code for extraction of coefficients is done in the while loop which exits when it encounter string "END" in the dat file.
4. Coming back to line1, strsplit command is used to split strings by space between them & names are stored in element_names.
5. In the same line1, to get the local temperatures, line1 is flipped using flip command & position 2, 3 & 4 are local mid , high , low temperatures respectively.
6. To get the coefficients, each coefficient are made up of 15 value using this information, A for loop is used to obatin this values & are stored in coef_mat.
7. Unnecessary variables are cleared whenever required.
8. For the plots , A range is required as we have two range one for low temperature & other for high temperature, A variable data_pts is assigned for that purpose.
9. The names of compounds are in string data type, but to use them in for loop they are converted into character data type.
10. The code for calculation of Cp, H & S is pretty straightforward ahead.
11. Lastly, the plots for all three properties for each compound is generated & saved in unique folder as required.

Main Code :

clear all;
close all;
clearvars;
clc;

% open a file descriptor
fd = fopen('THERMO.dat','r');

% Skip first line
fgetl(fd);

% Scan 2nd line temperature range
global_temp = (fscanf(fd, '  %f  %f  %f\n')).';

% Skip next 3 lines
fgetl(fd);
fgetl(fd);
fgetl(fd);

% Read each data point

element_names = [];
coef_mat = [];
local_temp = [];
line1 = fgetl(fd);
while ~strcmp(line1, 'END')
    % only read element name from 1st line
    a = strsplit(line1);
    element_names = [element_names; a(1, 1)];
    
    % for local temperatures
    b = flip(a);
    local_low = str2double(b{4});
    local_mid = str2double(b{2});
    local_high = str2double(b{3});
    local_temp = [local_temp; local_low, local_mid, local_high];
    % local_temp - end

    % extract coefficients
    line2 = fgetl(fd);
    line3 = fgetl(fd);
    line4 = fgetl(fd);                      

    num_start_idx = (1:15:61).'; % . is used to make it column vector
    [m, ~] = size(num_start_idx); % why not just m - tanu
    coef_row = zeros(1, 15);

    for i = 1:m
        idx = num_start_idx(i);
        val_line2 = str2double(line2(idx:idx+14));
        val_line3 = str2double(line3(idx:idx+14));
        val_line4 = str2double(line4(idx:idx+14));
        coef_row(1, i) = val_line2;
        coef_row(1, i + 5) = val_line3;
        coef_row(1, i + 10) = val_line4;
    end 
    
    coef_mat = [coef_mat; coef_row];
    line1 = fgetl(fd);
end

fclose(fd);

% drop last col - NAN , Removing 15th entry from each elements
coef_mat = coef_mat(:, 1:14);

% clear unwanted variables
% remove when using this as function
clear fd  line1 line2 line3 line4
clear i a idx m num_start_idx val_line2 val_line3 val_line4
clear ans
clear coef_row
clear b local_low local_mid local_high 




% Calculation of constants

% local_temp(:,1) = low temperature range
% local_temp(:,3) = high temperature range
% coef_mat column 1 to 7 are high temperature coefficients
% coef_mat column 8 to 14 are low temperature coefficients

R = 8.314472; % Declaration of Universal Gas Constant


 data_pts = 100; % It represents the number of data points between temperature range
 cp = zeros(53, data_pts*2);
 h = zeros(53, data_pts*2);
 s = zeros(53, data_pts*2);
 temperature = zeros(53,data_pts*2);
 
 species = char(element_names); % Converting string into character data type

 
 for j = 1:length(element_names)
     
     t_l = linspace(local_temp(j,1),local_temp(j,2),data_pts); % low temperature range
     t_h = linspace(local_temp(j,2)+1,local_temp(j,3),data_pts);%high temperature range
     
     % CP - Specific Heat

     cp_low = zeros(1,data_pts);
     cp_high = zeros(1,data_pts);


     for i=1:data_pts

        cp_low(1,i) = R.*( coef_mat(j,8) ...
                + coef_mat(j,9).*(t_l(1,i))...
                + coef_mat(j,10).*(t_l(1,i)).^2 ...
                + coef_mat(j,11).*(t_l(1,i)).^3 ...
                + coef_mat(j,12).*(t_l(1,i)).^4);

        cp_high(1,i) = R.*( coef_mat(j,1) ...
                + coef_mat(j,2).*(t_h(1,i))...
                + coef_mat(j,3).*(t_h(1,i)).^2 ...
                + coef_mat(j,4).*(t_h(1,i)).^3 ...
                + coef_mat(j,5).*(t_h(1,i)).^4);
     end

     cp(j,:) = [cp_low cp_high];
     
     
      % H - Enthalpy
     
      h_low = zeros(1,data_pts);
      h_high = zeros(1,data_pts);


     for i=1:data_pts

        h_low(1,i) = (R.*t_l(1,i)).*( coef_mat(j,8)...
                + (coef_mat(j,9).*(t_l(1,i)))./2 ...
                + (coef_mat(j,10).*(t_l(1,i)).^2)./3 ...
                + (coef_mat(j,11).*(t_l(1,i)).^3)./4 ...
                + (coef_mat(j,12).*(t_l(1,i)).^4)./5 ...
                + (coef_mat(j,13)./t_l(1,i)) ) ;
            
        h_high(1,i) = (R.*t_h(1,i)).*(coef_mat(j,1)...
                + (coef_mat(j,2).*(t_h(1,i)))./2 ...
                + (coef_mat(j,3).*(t_h(1,i)).^2)./3 ...
                + (coef_mat(j,4).*(t_h(1,i)).^3)./4 ...
                + (coef_mat(j,5).*(t_h(1,i)).^4)./5 ...
                + (coef_mat(j,6)./t_h(1,i)) ) ;

        
     end

     h(j,:) = [h_low h_high];
     
     % S - Entropy
     
     s_low = zeros(1,data_pts);
     s_high = zeros(1,data_pts);


     for i=1:data_pts

        s_low(1,i) = R.*( (coef_mat(j,8) .* log(t_l(1,i))) ...
                + (coef_mat(j,9).*(t_l(1,i)) ) ...
                + (coef_mat(j,10).*(t_l(1,i)).^2)./2 ...
                + (coef_mat(j,11).*(t_l(1,i)).^3)./3 ...
                + (coef_mat(j,12).*(t_l(1,i)).^4)./4 ...
                + coef_mat(j,14) ) ;

        s_high(1,i) = R.*( (coef_mat(j,1) .* log(t_h(1,i))) ...
                + (coef_mat(j,2).*(t_h(1,i)) ) ...
                + (coef_mat(j,3).*(t_h(1,i)).^2)./2 ...
                + (coef_mat(j,4).*(t_h(1,i)).^3)./3 ...
                + (coef_mat(j,5).*(t_h(1,i)).^4)./4 ...
                + coef_mat(j,7) ) ;
                
     end

      s(j,:) = [s_low s_high];
      
     
     temperature(j,:) = [t_l t_h]; 
     
     % Plot Generation
     
     % Creating Folder
    folder_name = sprintf( '%d)  %s',j,species(j,:)); 
    current_dir = pwd;
    mkdir(folder_name); % Create folder
    cd(folder_name); % To open that folder

    figure(1);
        plot(temperature(j,:),cp(j,:),'r','linewidth',3);
        grid on
        xlabel('Local Temperature (K)')
        ylabel('Specific Heat (J/kg-K)')
        title('Variation of Specific Heat w.r.t Temperature')
    image1 = sprintf("Cp Vs. T of %s.png",species(j,:)); % naming file
    saveas(gca,image1) % saving image

    figure(2);
        plot(temperature(j,:),h(j,:),'b','linewidth',3);
        grid on
        xlabel('Local Temperature (K)')
        ylabel('Enthalpy (J)')
        title('Variation of Enthalpy w.r.t Temperature')
    image2 = sprintf("H Vs. T of %s.png",species(j,:)); % naming file
    saveas(gca,image2) % saving image

    figure(3);
        plot(temperature(j,:),s(j,:),'g','linewidth',3);
        grid on
        xlabel('Local Temperature (K)')
        ylabel('Entropy (J/K)')
        title('Variation of Entropy w.r.t Temperature')
    image3 = sprintf("S Vs. T of %s.png",species(j,:)); % naming file
    saveas(gca,image3) % saving image

    cd(current_dir); % Returning to current directory
     
     % Plot - end
     
 end
 % Clearing unwanted variables to mantain workplace minimal
 clear cp_low cp_high h_low h_high i j s_low s_high t_l t_h

Code for calculation of Molecular Weight :

1. For the purpose of calculation of molecular weight, MATLAB (containers)mapping is used.
2. In mapping a key_element is assigned a value, for eg here 'C' for carbon atom its molecular weight is assigned '12.0107'.
3. Next, The for loop is initiated to run for all 53 compunds/elements.
4. Using strip() command the empty spaces are eliminated.
5. So that the while loop runs only for the number of characters present in that string.
6. Each character is compard with the mapped elements & their respective molecular weight is appened, When a digit is encountered it is stored alongwith its previous character.
7. At the end the final molecular weight is calculated, stored & displayed in command window.

element = {'H','C','N','O','A','(','S',')','R'};
atomic_weight = [1.0078 12.0107 14.0067 15.9949 39.9480, 0.00, 32.0600, 0.00, 0.00];
elem_weight_map = containers.Map(element, atomic_weight);

molecular_weight = zeros(length(element_names),1);
for jj = 1:length(element_names)
    
    speci = strip(species(jj,:));
    %fprintf('%s\n', speci);

    idx = 1;
    len = length(speci);
    
    while idx <= len
    
        elem = speci(1, idx);
   
        num = 1;
        if (idx < len) && isstrprop(speci(1, idx+1), 'digit')
            num = str2double(speci(1, idx+1));
            idx = idx + 2;
        else
            idx = idx + 1;
        end
    
        molecular_weight(jj,1) = molecular_weight(jj,1) + elem_weight_map(elem)*num;

    end
    fprintf('%2.0f.The molecular weight of %s is %f. \n \n ',jj,speci,molecular_weight(jj,1))
end

Molecular Weights (In amu) :

Plots of all species :

Plots for (`O_2`)

Plots for (`N_2`)

Plots for (`CO_2`)

Erros Faced : Several errors were faced , but the major ones were the following:

1. Incorrect parenthesis in the calculation of Cp, H & S.

2. To get the value of local temperatures firstly I find the postion of G & next three values would be the temperature values, But for some elements this caused a problem because the number of elements in that row were different from the rest, Therefore to avoid that another logic was used incorporationg flip().

Conclusion : During the B.Tech course generally, The value of above properties are taken constant for various problem solving which does not hold true for real-world application, Hence the incorporation of such variation are important.

With the help of Data parsing, these variation of properties can be esaily used in mathematical problem solving in MATLAB.