Parsing NASA Thermodynamic Data

Parsing NASA Thermodynamic Data

1. AIM:

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

2. GOVERNING EQUATIONS IF ANY:

1. `C_p = R*(a_1 + a_2T + a_3T^2 + a_4T^3 + a_5T^4)`

2. `H = R*T*(a_1 + a_2T/2 + a_3T^2/3 + a_4T^3/4 + a_5T^4/5 + a_6/T)`

3. `S = R*(a_1*ln(T) + a_2T + a_3T^2/2 + a_4T^3/3 + a_5T^4/4 + a_7)`

Here, R = Universal Gas Constant (`J/(mol*K)`)

H = Enthalpy of a specie (`(KJ)/(Kg)`)

S = Entropy of a specie (`J/K`)

3. OBJECTIVES OF THE PROJECT:

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

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

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

4. FLOW OF THE PROJECT:

1. Opening the data file and check for key things like unnecessary lines and line of start i.e from where actual data can be extracted.

2. Extracting global and local minimum, medium and maximum temperatures.

3. Looping to extract all the 14 coefficients line wise and assigning then to a variable.

4. Using governing equations, finding the values of Specific heat, Enthalpy and Entropy.

5. Plotting the graphs for `O_2 , N_2 and CO_2`.

6. Saving all plots in individual folder.

7. Calculating molecular weight and displaying it in command window.

5. MATLAB CODE:

Main program (Only THERMO.dat as input and plots, weights as outputs)

clear all
close all
clc

f1 = fopen('THERMO.dat','r');            %open THERMO.dat file with read only permission
fgetl(f1);                               %print the first line
b = str2num(fgetl(f1));                  %convert the string in the next line into numbers
g_low_temp = b(1);                       %assign values of numbers from previous statement to variable b()
g_med_temp = b(2);
g_high_temp = b(3);
as = 'END';
r = 8.314;                               %Universal gas Constant J.mol^-1 K^-1

output = main_func(f1,as,r);             %Function to get output of Cp, H and S

Code Explaination:

- We have opened the file named THERMO.dat giving only read permission.

- We have skipped lines so that we can extract global temperatures i.e. Low, Medium and High from the second line and assigned it to variables using string to number conversion.

- Created a variable named as = 'END' so as to determine the end of the dat file i.e. how many lines are there in the dat file. That will be used later.

main_func(f1,as,r)

function [final,weight] = main_func(f1,as,r)

for i=1:500                              %For loop to find how many lines are there in the THERMO.dat file
    a = fgetl(f1);                       %Advance to next line with every iteration
    if strcmp(a,as)==1                   %For loop to compare as = 'END' with the last line of THERMO.dat file
        c = i+2;
        break
    end
end
    
    for j=1:c                           %For loop to come back to the line after which we started previous for loop
        frewind(f1);                    %Rewind line function so that we can start from initial lines
    end
    
loop_count = (i+2-5-1)/4;               %Determining how many species are present. -5 and -1 are lines which are not needed.

for i=1:5                               %Skipping unwanted lines
    fgetl(f1);
end

for start=1:loop_count                  %For loop to run the code for the as many times as the number of species
c = fgetl(f1);
species = strtok(c)                     %Finding the species by strtok i.e all strings before the default delimiter ' '
d = findstr(c,'G ');                    %Every local temperature is preceded by this string 'G '. Thus find it.
f = c(d+1:length(c));                   
e = str2num(f);                         %Convert the local temp string f into number
l_low_temp = e(1);                      %Assign the temp to variable e
l_high_temp = e(2);
l_med_temp = e(3);

q = fgetl(f1);
g = findstr(q,'E');                     %Checkpoint to distinguish 14 coefficients from a single line
int = 1;

for i=1:5                               %Since 5 coefficients are there in a single line, i=5 . Starting from line number 2
    h = q(int:g(i)+3);                  %Take coefficient from initial position to position of E + 3 positions.
    aa(i) = str2num(h);
    int = g(i)+4;                       %take the cursor before next coefficient identification
    if q(int) == ' '                    %If the coefficient is not negative, there will be a space before it. 
        int = int+1;                    %Thus IF condition so that we can skip that space.
    else
        continue
    end
end

q = fgetl(f1);
g = findstr(q,'E');
int = 1; 

for i=1:5                               %for loop for coefficient identification for line number 3
    h = q(int:g(i)+3);
    ab(i) = str2num(h);
    int = g(i)+4;
    if q(int) == ' '
        int = int+1;
    else
        continue
    end
end

q = fgetl(f1);
g = findstr(q,'E');
int = 1;

for i=1:4                                %for loop for coefficient identification for line number 4
    h = q(int:g(i)+3);
    ac(i) = str2num(h);
    int = g(i)+4;
    if q(int) == ' '
        int = int+1;
    else
        continue
    end
end

final = specific(aa,ab,ac,e,r,species);  %Function output for plots and writing image files
weight(start) = molweight(species)      %Function to display molecular weights in notepad

end
end

Code Explaination:

- Here we first have created a loop in order to determine how many lines are there in the dat file. For finding out number of lines, we have used fgetl() function en-number of times. In order to restore the location of fgetl(), we have used the frewind() function so that line cursor can be bought to original location.

- loop count is nothing but the number of species present in the dat file. +2 is the first 2 lines which we have already skipped using fgetl(), -5 is the number of header lines which we can skip, -1 is the last line which we can skip and divided by 4 ( /4 ) is the number of lines present in each specie. Thus total species is obtained so that we can loop that many number of times.

- Skip 5 lines to the start of specie #1

- Using strtok, we can extract all characters of string till we reach a delimiter which is ' ' by default. Thus we get the name of the specie.

- Then we use findstr() to locate a checkpoint, in this case it is 'G ' so that we can extract local temperatures which are present after 'G '.

- We use str2num() to convert the string into numbers and then assign them to variable e()

- In order to extract the 14 coefficients, 5 present on each line, we need such a checkpoint which is 'E' present in every coefficient.

- Thus using findstr() we find position / occurance of 'E' in one line and then according to that, we take all the preceding numbers and 3 numbers after 'E' so that we get a complete coefficient.

- Same logic is followed for next 2 lines and it is also explained in the code using comments.

- If there is no negative sign, the next coefficient is preceded by a space ' '. That is also taken care of by using IF condition.

- At the end, 2 more functions are used. They are as follows:

1. specific() is used to plot the graphs and save them

2. molweight() is used to calculate molecular weight and display them in command window.

specific(aa,ab,ac,e,r,species)

function [cp1,cp2,h1,h2,s1,s2] = specific(aa,ab,ac,e,r,species)

x = linspace(e(1),e(3),200);              %200 temperature values between local minimum and local medium temps
y = linspace(e(3),e(2),200);              %200 temperature values between local medium and local maximum temps

cp1 = r*(aa(1)+ x.*aa(2)+ x.^2.*aa(3) + x.^3.*aa(4) + x.^4.*aa(5));
cp2 = r*(ab(3)+ y.*ab(4)+ y.^2.*ab(5) + y.^3.*ac(1) + y.^4.*ac(2));
figure(1)
plot(x,cp1,'linewidth',2);
hold on
plot(y,cp2,'linewidth',2);
hold off
xlabel('Temperature')
ylabel('Specific Heat')

x1 = 'D:\Folders\SKILL-LYNC\Matlab\Challenges\nasa_therm_data\plots';    %Navigate to this directory
cd(x1);                                  %change directory to above mentioned directory
mkdir(species);                          %create directory named after the respective species
cd(species);                             %change directory to above created directory
saveas(gcf,'Specific Heat')              %Save the plot as .fig
cd(x1);                                  %again go back to parent directory x1.

h1 = x.*(aa(1)+ x.*aa(2)./2+ x.^2.*aa(3)./3 + x.^3.*aa(4)./4 + x.^4.*aa(5)./5 + ab(1)./x)*r;
h2 = y.*(ab(3)+ y.*ab(4)./2+ y.^2.*ab(5)./3 + y.^3.*ac(1)./4 + y.^4.*ac(2)./5 + ac(3)./y)*r;
figure(2)
plot(x,h1,'linewidth',2);
hold on
plot(y,h2,'linewidth',2);
hold off
xlabel('Temperature')
ylabel('Enthalpy')

cd(species);                             %No need to create directory again as it is created earlier.
saveas(gcf,'Enthalpy')
cd(x1);

s1 = (aa(1)*(log(x)/log(exp(1))) + x.*aa(2) + x.^2.*aa(3)./2 + x.^3.*aa(4)./3 + x.^4.*aa(5)./4 + ab(2))*r;
s2 = (ab(3)*(log(y)/log(exp(1))) + y.*ab(4) + y.^2.*ab(5)./2 + y.^3.*ac(1)./3 + y.^4.*ac(2)./4 + ac(4))*r;
figure(3)
plot(x,s1,'linewidth',2);
hold on
plot(y,s2,'linewidth',2);
hold off
xlabel('Temperature')
ylabel('Entropy')

cd(species);
saveas(gcf,'Entropy')
cd(x1);
end

Code Explaination:

- The local temperatures are divided into 2 temperature ranges by linspace(). One range is from Low to Medium and second range is from Medium to High.

- The function then contains formulas of governing equations.

- Next, we open a new figure and plot the graph of local temperature and Cp.

- Next we navigate to a directory where my code is stored. Then we change the directory to 'plots' which is inside my parent directory. Then we make a new directory whose name is the same as the specie being analyzed. and save the plot inside that newly made directory. Then we again change the directory to parent directory.

- Same steps are repeated for Enthalpy as well as Entropy.

molweight(species)

function [mol] = molweight(species)

    aname = ['H','C','O','N','AR'];                      %string of parent elements
    aweight= [1.008 12.011 15.999 14.007 39.948];       %array of molecular weights of respective parent elements
    moll = 1;
    for i=1:length(species)                             %To run the loop as many times as length of the species                    
  
        if str2num(species(i))>1                        %to identify if the object at species(i) is string or a number
            moll = moll + (str2num(species(i))-1)*aweight(j);
            %If it is a number, just multiply it by str2num -1 because
            %already molecular weight is added in previous string
        else
            j=1;
            while strcmp(species(i),aname(j)) ~= 1 && j<5   %to check which parent element is present
                 j=j+1;
            end
            moll=moll+aweight(j);
            
        end
    end
    
    mol = moll-1;                %subtracting initial moll = 1
    
    f1 = fopen('aaa.txt','a+');  %creating a new text file
    fprintf(f1,'Molecular weight of %s is %f \n',species,mol)     %writing molecular weight of element on new line
    fclose(f1);                  %closing the file. the attribute a+ is for read and write while keeping previous content intact.
end

    

Code Explaination:

- Since every compound is made of some parent element, we have pre-defined some parent elements H, C, O, N, Ar along with their molecular weights.

- Then in order to determine if every string letter of a specie is an alphabet or a number, we have used str2num() function.

- We do not write 1 infront of any element. For eg: We do not write `N_1H_3`. Instead we write `NH_3`. So we can say that if the str2num() of any string letter is more than 1, then it is definitely a number. Else it is an alphabet.

- The for loop then adds the molecular weight accordingly.

- If it isn't a number then the while() loop finds the matching element and then adds the molecular weight likewise.

- At the end, we have created a text file having attributes Read and write without deleting existing content and then we have printed the molecular weight of each specie in it linewise.

6. OBSERVATIONS:

For `N_2`,

For `O_2`,

For `CO_2`,

The folder dumps for each and every specie:

The Notepad containing molecular weights:

The command window having all the molecular weights:

7. ERRORS FACED:

Fixed by using strcmp() so that string can be compared.

- Error where the plot could not be saved as a png or jpeg or any format. The solution to this error was not found after many attempts and thus, the plots had to be saved as a '.fig' format which is a matlab plot format.

- Rest errors were of invalid syntax and those were corrected on the spot.

8. CONCLUSION:

Thus we have successfully parsed the thermodynamic data from a data file and plotted the graphs of different thermodynamic properties and calculated the molecular weights of different species present in the data file.