Parsing NASA thermodynamic data

1. AIM

• NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S. They have also documented the co-efficients that are required to evaluate these polynomials for 1000+ species.

2. OBJECTIVE

• 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.
• Calculate the molecular weight of each species and display it in the command window.
• Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species.
• Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species.

3. BODY OF THE CONTENT

• Definition of File parsing:- The definition of parsing is “Analyzing the data”, in this case, we are analyzing the data from “THERMO.DAT” file, and calculating useful information.

Logic Behind Main program coding Coding: -

Formulae Used;-

• In this case, temperature is varying, hence for first 7 coefficient we use local high temperature & for rest of coefficient we will use local low temperature.
• To calculate molecular weight we will use

Molucular weight=(Number of items)*(Molecular weight)

4. Program Explanation:-

Main program(“File_Parsing”):-

close all;
clear all;
clc;

%Opening the file
file = fopen('THERMO','r');  

%Getting information of species from function
 for i=1:53;

     [Cp,H,S,Temperature,species_name]=coefficients(file,i);   %Calling the function to calculate Cp,H & S
     Molecular_mass = molecular_mass(species_name);          %Calling function to calculate Molucular_mass
     disp(['Molucular weight of ',species_name])
     disp(Molecular_mass)

%plotting
%Temperature vs Entropy
figure(1)
plot(Temperature,S,'linewidth',4,'color','b')
title(species_name)
xlabel('Temperature')
ylabel('Entropy')

%Temperature vs Enthalpy
figure(2)
plot(Temperature,H,'linewidth',4,'color','k')
title(species_name)
xlabel('Temperature')
ylabel('Enthalpy')

%Temperature vs Specific heat
figure(3)
plot(Temperature,Cp,'linewidth',4,'color','r')
title(species_name)
xlabel('Temperature')
ylabel('Specific heat')

%changing the filepath to save the plots
cd(species_name)

%saving the plots

saveas(1,'Entropy.jpg');
saveas(2,'Enthalpy.jpg');
saveas(3,'Specific heat.jpg');

cd('C:\study\');
 end

• In this main program, as explained earlier we are calling the functions for calculating value of Cp, H and S. Then we plot the graphs for each species.
• In line no. 6, we are opening the file “THERMO”, in which all data is stored.
• From line no. 9, we are using for loop, as here are total 53 species, and we have to repeat same steps of calculating coefficients & Molecular weight.
• In line no. 11 we are calling function for calculating Cp, H and S, in subsequent line we are calling function for calculating molecular weight.
• In line no. 13 & 14 we use “disp” command to show molecular weight in specific format.
• From line no. 17 to 36, we are plotting 3 graphs for Temperature vs Entropy ,Temperature vs Enthalpy, Temperature vs Specific heat respectively.
• In line no. 39 we are creating separate folder with “species name”, to save graph in it using “saveas”command,(line no. 43 to 45)
• In line no 47 we are changing bach directory to its precious location to continue the loop againg& ending the loop with “end ” command.

Program for calculating Cp, H and S:-

function [Cp,H,S,Temperature,species_name,Molecular_mass]=coefficients(file,i)

if i==1;
    line1=fgetl(file);
line2=fgetl(file);

% Ignorring next three line, as they have no use in programming point of
% view
comment_1=fgetl(file);
comment_2=fgetl(file);
comment_3=fgetl(file);

X=strsplit(line2,' ');
Globle_lower_temp=str2num(X{2});
Globle_middle_temp=str2num(X{3});
Globle_high_temp=str2num(X{4});
end

     species=fgetl(file);
     species_name=species(1:6);
     l1=fgetl(file);
     l2=fgetl(file);
     l3=fgetl(file);
     
     %Local temperature rage
     t = findstr(species,'.'); %To find the location of "."

% Local Temperature for each species
t_low=species(52-3:52+3);
Local_Temperature_low=str2num(t_low);

t_mid=species(72-4:72+3);
Local_Temperature_mid=str2num(t_mid);

t_high=species(62-4:62+3);
Local_Temperature_high=str2num(t_high);


Temperature =linspace(Local_Temperature_low,Local_Temperature_high,100);


     %Assigning coefficients
     %Finding position of E, to extrack coefficients.
     line_1=strfind(l1,'E');
     line_2=strfind(l2,'E');
     line_3=strfind(l3,'E');
     
     c1=str2num(l1(1:line_1(1)+3));
     c2=str2num(l1((line_1(1)+4):(line_1(2)+3)));
     c3=str2num(l1((line_1(2)+4):(line_1(3)+3)));
     c4=str2num(l1((line_1(3)+4):(line_1(4)+3)));
     c5=str2num(l1((line_1(4)+4):(line_1(5)+3)));
     
     c6=str2num(l2(1:line_2(1)+3));
     c7=str2num(l2((line_2(1)+4):(line_2(2)+3)));
     c8=str2num(l2((line_2(2)+4):(line_2(3)+3)));
     c9=str2num(l2((line_2(3)+4):(line_2(4)+3)));
     c10=str2num(l2((line_2(4)+4):(line_2(5)+3)));
     
     c11=str2num(l3(1:line_3(1)+3));
     c12=str2num(l3((line_3(1)+4):(line_3(2)+3)));
     c13=str2num(l3((line_3(2)+4):(line_3(3)+3)));
     c14=str2num(l3((line_3(3)+4):(line_3(4)+3)));
     
     %creating folder for each species
    mkdir(species_name)

     % Caluclulation of Cp, Entropy & Enthalphy
     
     R =8.314;                              %Universal gas constant

     %calculation of properties based on the lower and upper temperature

    if Temperature<3000
Cp = (c1+c2.*Temperature+c3.*Temperature.^2+c4.*Temperature.^3+c5.*Temperature.^4).*R;
H = (c1+(c2.*Temperature)./2+(c3.*Temperature.^2)/3+(c4.*Temperature.^3)/4+(c5.*Temperature.^4/5)+c6./Temperature).*Temperature.*R;
S = (c1*log(Temperature)+c2.*Temperature+(c3.*Temperature.^2)/2+(c4.*Temperature.^3)/3+(c5.*Temperature.^4)/4+c7).*R;
    else
Cp = (c8+c9.*Temperature+c10.*Temperature.^2+c11.*Temperature.^3+c12.*Temperature.^4).*R;
H = (c8+(c9.*Temperature)/2+(c10.*Temperature.^2)/3+(c11.*Temperature.^3)/4+(c12.*Temperature.^4/5)+c13./Temperature).*Temperature.*R;
S = (c8.*log(Temperature)+c9.*Temperature+(c10.*Temperature.^2)/2+(c11.*Temperature.^3)/3+(c12.*Temperature.^4)/4+c14).*R;
  
end

• In line no. 11 of main program we are calling the function “coefficient” for calculating Cp, H and S.
• In this program, we are analising previously opened file “THERMO”.
• From line no. 3 to 17 we are creating loop, as we needed information like file name, global temperature range only at first time, furthermore there is no use of that information.
• As we know, the contain of file is in “string” format, hence we cannot analize it unless we convert it to numeric format.
• Hence we we use “strsplit” & “str2num” command.
• We know at line no. 6 of “THERMO” file, there are global temperature rang of species is given, and those are separated by “ (space)”.
• Hence first we use “strsplit” command which devide the string in 3 srings by identifying spece between them.
• Then we convert these three individual strings to numeric by using “str2num” for calculation purpose.
• In line no. we are identifying respective species name for calculation.
• As we have to use local temperature not globle temperature for caluculation, we have to file local temperature range which is given in same line in original file as in species name is given.
• Local temperature is separated by “.”, hence we use “findstr” to file location of “.” In given line, then we use “str2num” to convert given string to numeric format.
• By this technique from line no. 29 to 39. We find local low, mid & high temperature. & put it together in matrix format for plotting the graph.
• In subsequent lines we have given coefficents in “THERMO” file, but to extraxt them we have to build logic, for this we know that each coeffient contains letter “E”, after which there are three letters then next coefficient starts.
• Hence first we fine location of letter “E” in line using “strfind”, then we convert this string to numeric.
• By using this logic we fing 14 coefficent from line no 42 to 63.
• Then we find the value of Cp, H and S, by using formula as stated previously(line no. 73 to 82)
• We used loop here calucalting no. of values of Cp, H and S for each temperature, also we use “if” command for logic that, if temperature is below 3000 then use coefficient from c1 to c7 else used coefficient from c8 to c14.

Program for calculating Molucular weight:-

function Molucular_weight = molecular_weight(species_name)
  
species = upper(species_name); 

Molucule = ['H', 'C', 'N', 'O', 'A'];                   % Making consttituents name

Atomic_weight = [1.00794 12.0107 14.0067 15.9994 39.948];   %Creating array for molucular weight

Molucular_weight = 0; 

for i=1:length(species)                                     %Creating loop for comparring
    for j=1:length(Molucule)
        if strcmp(species(i),Molucule(j))               %Comparing species name with constituents
           Molucular_weight = Molucular_weight + Atomic_weight(j);%if both are same then calculatingmolucular weight
           u = j;
        end
    end
    n = str2num(species(i));                
    if n>1                                                  %if species contain more then one moluculem, then recalculate molucular weight
        Molucular_weight = Molucular_weight + Atomic_weight(u)*(n-1);
    end
end
end

• We are taking species name as input & molecular weight as output.
• We made two matrix, molecule & Atomic_weight, which contains total molecules & respective molecular weight.
• Initially we set molecular weight to 0, then we compare input species with each element of “molecule” matrix. If it match then we take respective molecule weight of respective molecule & calculate molecular weight of species. (line 9 to 17)
• If no. of molecule is more than 1 in given species then we add previous value of molecular weight & recalculate it. (line no. 18)

5.Result

Folder for all save plots:-

Command windows which showing molecular weight of species: -

Plot of O2:-

Plot for N2:-

Plot for Co2:-

Errors faced during program:-

Straight line graph

Remedy:- previously I was using, only two values of temperature(t_low & t_high) ,hence getting only two values, so graph is straight line. now I used “linpsace” command to use 100 values of temperature between given lower temperature to local high temperature.