Project 1 - Parsing NASA thermodynamic data

Aim:

To write a function that extracts the 14 coefficients 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.

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

Formula used:

`C p = R ( C 1 + C 2 ⋅ T + C 3 ⋅ T 2 + C 4 ⋅ T 3 + C 5 ⋅ T 4 )` ------> Specific heat

`H = R T ( C 1 + C 2 ⋅ T 2 + C 3 ⋅ T 2 3 + C 4 ⋅ T 3 4 + C 5 ⋅ T 4 5 + C 6 T )` ---------> Enthalpy

`S = R ( C 1 ⋅ ln ( T ) + C 2 ⋅ T + C 3 ⋅ T 2 2 + C 4 ⋅ T 3 3 + C 5 ⋅ T 4 4 + C 7 )` ----------> Entropy

Here R = 8.314 is the universal gas constant, T is the temperature.

Main code:

clear all
close all
clc

%Inputs
R= 8.314;

%Reading the data file
file = fopen('THERMo.dat','r');

%Reading the first line of the file 
l1 =fgetl(file);

%Getting the global minimum, maximum and mid temperatures
l2=fgetl(file);
temperatures = strsplit(l2);
global_min = str2num(temperatures{2});
global_mid = str2num(temperatures{3});
global_max = str2num(temperatures{4});

%Defining the temperature range 
temp = linspace (global_min,global_max,500);

%Ignoring the comments in the file 
for k=1:3
fgetl(file);
end

%Creating a loop to go through all the 53 species of data in the file 
for j=1:53
    %Reading the first line which includes the information for species
    l=fgetl(file);
    
    %Reading the name of the species
    element_name=strtok(l);
    
    %co_efficients
    l1=fgetl(file);
    l2=fgetl(file);
    l3=fgetl(file);
    
    %Finding the position of E so that we can extract the coefficients
    a=strfind(l1,'E');
    b=strfind(l2,'E');
    c=strfind(l3,'E');

    % high temp coeff
    T1=str2num(l1(1:a(1)+3));
    T2=str2num(l1(a(1)+3:a(2)+3));
    T3=str2num(l1(a(2)+4:a(3)+3));
    T4=str2num(l1(a(3)+4:a(4)+3));
    T5=str2num(l1(a(4)+4:a(5)+3));
    T6=str2num(l2(1:b(1)+3));
    T7=str2num(l2(b(1)+4:b(2)+3));

    %Extracting low temp coeff
    T8=str2num(l2(b(2)+4:b(3)+3));
    T9=str2num(l2(b(3)+4:b(4)+3));
    T10=str2num(l2(b(4)+4:b(5)+3));
    T11=str2num(l3(1:c(1)+3));
    T12=str2num(l3(c(1)+4:c(2)+3));
    T13=str2num(l3(c(2)+4:c(3)+3));
    T14=str2num(l3(c(3)+4:c(4)+3));

    %Calculating specific heat
    c_p = specificheat(T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,R,temp,global_min,global_max,global_mid);
    
    %Calculalting enthalpy
    H = enthalpy(T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,R,temp,global_min,global_max,global_mid);
    
    %calculating entropy
    S = entropy(T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,R,temp,global_min,global_max,global_mid);
    
    %Plotting graphs
    plotting(c_p,H,S,temp,element_name)
    
    %Calculating molecular weights
    mol_wt(j) = molecularmass(element_name)
    
    %Function to write molecular weight to a file 
    weight = mol_wt(j);
    write_mol_weight(element_name, weight)
    
end

fclose (file)

• The value of the universal gas constant R is given as 8.314 as an input.
• The "fopen" command is used to open. Then the file is read by specifying 'r' keyword in the "fopen" command.
• The file is read line by line using the "fgetl" command. The first line represents the title so it is ignored.
• The next line has the global maximum, global mid, and global minimum values which are read by using "fgetl" command and split into the three values by using the "strsplit" command and converted into numbers.
• Then the temperature array is defined using "linspace" command as 500 values between the global minimum and global maximum values.
• Then "for" loop is created to loop through the 53 species data given in the file and calculate the required values.
• Inside the loop the "fgetl" command is used to read the next line of the file which represents the name of the species.
• The name of the element is stored in a variable using strtok command.
• Then the next 3 lines of the file are read which contains the high and low-temperature data of the particular species.
• The read lines are in the form of a string so they are converted into numbers using the str2num command and stored in specific variables so they can be used for further calculations.
• The low temperatures values are stored in variables T1, T2, T3, T4, T5, T6, and T7, and the high-temperature values are stored in variables T8, T9, T10, T11, T12, T13, T14.
• The next line is the reference to a function called specific heat which calculates and returns the value of specific heat using the high and low-temperature values as inputs.
• The next line is the reference to a function called enthalpy which calculates and returns the value of enthalpy using the high and low-temperature values as inputs.
• The next line is the reference to a function called entropy which calculates and returns the value of entropy using the high and low-temperature values as inputs.
• The next line is the reference to a function that plots graphs using the given inputs.
• The next line is the reference to a function that calculates the molecular weights of the different elements.
• The next line is the reference to a function that writes the molecular weights of different elements to a file.
• The NASA thermodynamic data file is closed using the "fclose" command.

code to calculate Specific heat:

function c_p = specificheat(T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,R,temp,global_min,global_max,global_mid)
    if temp<=global_mid
        c_p = R*(T1 + T2.*temp + T3.*(temp.^2) + T4.*(temp.^3) + T5.*(temp.^4))
    else
        c_p = R*(T8 + T9.*temp+T10.*(temp.^2) + T11.*(temp.^3) + T12.*(temp.^4))
end

we use the "if-else" conditional statement to check if the temp values are less than or equal to the global mid-temperature values and use appropriate equations with high and low-temperature coefficients to calculate the specific heat value denoted by c_p and return that value.

code to calculate enthalpy:

function H = enthalpy(T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,R,temp,global_min,global_max,global_mid)
    if temp <=global_mid
        H = R.*temp.*(T1 + T2.*(temp./2) + T3.*((temp.^2)./3) + T4.*((temp.^3)./4) + T5.*((temp.^4)./5) + T6./temp)
    else
        H = R.*temp.*(T8 + T9.*(temp./2) + T10.*((temp.^2)./3) + T11.*((temp.^3)./4) + T12.*((temp.^4)./5) + b6./temp)
end

By using the "if-else" conditional statement the temp values are compared with the global mid-temperature values and the appropriate equation with the proper maximum or minimum temperature coefficients is selected and the enthalpy is calculated and its value represented by H is returned.

code to calculate entropy:

function S =entropy(T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,R,temp,global_min,global_max,global_mid)
    if temp<=global_mid
        S = R*(T1.*log(temp) + T2.*temp + T3.*(temp.^2)./2 + T4.*(temp.^3)./3 + T5.*(temp.^4)./4 + T7)
    else
        S = R*(T8.*log(temp) + T9.*temp + T10.*(temp.^2)./2 + T11.*(temp.^3)./3 + T12.*(temp.^4)./4 + T14)

end

Then using the "if-else" conditional statement the temp values are compared with the global mid-temperature values and the appropriate equation with the proper maximum or minimum temperature coefficients is selected and the entropy is calculated and its value represented by S is returned.

code for plotting:

function plotting(c_p,H,S,temp,element_name)
cd(pwd)
dir_name=element_name;
mkdir(dir_name);

%Plotting temperature vs specific heat
figure(1)
plot(temp,c_p,'lineWidth',3)
title(sprintf('Specific Heat vs Temperature rn for %s:',element_name))
xlabel('Temperature')
ylabel('Specific Heat')

%Plotting temperature vs specific heat
figure(2)
plot(temp,H,'linewidth',3)
title(sprintf('Enthalpy vs Temperature rn for %s:',element_name))
xlabel('Temperature')
ylabel('Enthalpy')

figure(3)
plot(temp,S,'linewidth',3)
title(sprintf('Entropy vs Temperature rn for %s:',element_name))
xlabel('Temperature')
ylabel('Entropy')

cd(dir_name)
saveas(1,'c_p_vs_temp.png')
saveas(2,'H_vs_temp.png')
saveas(3,'S_vs_temp.png')

cd ..

end

Using the "mkdir "command a new folder is created in the current folder "pwd" with the same name as the element name.

The first figure(1) represents the plots the temperature vs specific heat values.

The first figure(2) represents the plots of the temperature vs enthalpy values.

The first figure(3) represents the plots of the temperature vs entropy values.

Using the "saveas" command the three figures are saved in the new folder which is created by the element name.

Using the cd.. command the folder is changed back to the original folder.

The function runs for 53 times through the "for" loop and it creates separate folders for each species and saves all the plots in the folders.

code to calculate molecular mass:

function [w] = molecularmass(element_name)

s = {'H' 'C' 'N' 'O' 'AR'}; 
%atomic weight of each elements
atomic_wt = [1 12 14 16 40];                                                      
w = 0;                                                                       

for i = 1:length(element_name)  
 for j = 1:length(s)
    if strcmp(element_name(i), s(j))                           
    w= w + atomic_wt(j); 
    %storing the position of matched element
    position = j;                                                         
  end 
end 
   
n = str2num(element_name(i));
if n > 1
                                                      
w = w + atomic_wt(position)*(n-1);                                        
end 
end

A string array is created which contains the name of the atoms.

A variable atomic_wt is created which contains the atomic weight of respective atoms in the string array.

Then used "strcmp" command to compare the name of the element with the string array.

If there is the repetition of any element, "str2num" is used to find that and (n-1) is multiplied to atomic_wt of the element to calculate the molecular weight.

code to write molecular mass:

    function write_mol_weight(element_name,weight)
    dir_name = element_name
    mol_file = fopen('mol_wt.docx' , 'a+')
    fprintf (mol_file, 'Molecular weight of %s is %f n', dir_name,weight)
    fclose (mol_file)
    end
 

The new file is created with the name 'mol_wt.docx' and it is opened in append mode denoted by 'a+' so that we can write the molecular weight of all elements in a single file.

Using the "fprintf" command the name of the element along with their molecular weight is written to the file.

The file is then closed.

The link for the file is,

link---> https://drive.google.com/file/d/19Nw3U74IRC5B9PSKuyTwfDohOk_JPOBp/view?usp=sharing

Results:

Output for one element is shown below:

Specific heat:

Enthalpy:

Entropy:

image showing the element-wise results saved in separate folders: