Project 1 - Parsing NASA thermodynamic data

Aim : -

• 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 then save the plots as images and dump them in separate folders for each species wih suitable name.

Given

File link

NASA thermodynamic data

Formulae

`(Cp)/R = a1+a2*T + a3*T^2 + a4*T^3 + a5*T^4`

`H/(R*T) = a1 + (a2*T)/2 + (a3*T^2)/3 + (a4*T^3)/4 + (a5*T^4)/5 + (a6)/T`

`S/R = a1*lnT + a2*T+(a3*T^2)/2 + (a4*T^3)/3 + (a5*T^4)/4 + a7`

where

Cp = Specific heat at constant pressure

H = Enthalpy

S = Entropy

R = Universal gas constant [8.314 J / (mol·K)]

T = Temperature

Theroy

Cp = Specific heat at constant pressure represents the heat supplied to a unit mass of the systemfhet to raise its temperature through 1K, keeping the pressure constant.

Enthalpy (H) = A thermodynamic quantity equivalent to the total heat content of a system. It is equal to the internal energy of the system plus the product of pressure and volume.

Entropy (S) = A thermodynamic quantity representing the unavailability of a system's thermal energy for conversion into mechanical work, often interpreted as the degree of disorder or randomness in the system

Steps:-

1. Start with clear all ,close all and clc for clearing by previous data.
2. Open THERMO.dat file for reading with fopen command and read permission and save it to f1, f1=fopen('THERMO.dat','r')
3. If f1 is less than zero there is error in opening the file for this make a if statements to print error
4. Make a folder "Plots" with mkdir command to save all the plots in single folder
5. Read first line with fgetl .
   • fgetl read one line at a time for reading next line use fgetl agian
6. Get second line of global temperature using fgetl command again and save the line to global_temp
7. Split global temperature line according to spaces with a built in function strsplit and save it to a cell A
8. The second , third and forth elements of the cell A are global minimum , medium and maximum temperature written in string class convert it double class and then save it as global_min_temp, global_mid_temp and global_max_temp
9. Use fgetl three times to skip next three line of comments
10. Create a for loop for all 53 species of the file
    1. Use fgetl to read first line of the species witch contain species name and local temperature for that species
    2. Split line according to spaces with strsplit function and form a cell L1
    3. Get the first element of cell L1 witch is species name and save it as species
    4. The forth, third and the second from the last are local minimum , maximum and midium temperature written in string class convert it to double using str2num and then save it as local_min_temp , local_max_temp and local_mid_temp .
    5. Find the values of first five coefficient with the help of position of E and save it as a1,a2,a3,a4 and a5
    6. Similarly find the next 5 cofficient from line 3 and last 4 coefficient from last line
    7. Define a temperature range from local min temperature to local max temperature containing 500 values and save it to a variable 'T' using linspace
    8. Make function thermodynamic_properties for calculation of Cp , H and S .
       1. Take all 14 coefficient , local_mid_temp and T as input to the function
       2. Define a universal gas constant R equal to 8.314 .
       3. Make a for loop with i goes from 1 to length of T .
          1. In the for loop if temperature is greater than local_mid_temp use first 7 coefficient for calculation of Cp , H and S using above given equation
          2. else if temperature is less tham local_mid_temp use last 7 coefficienet for calculation
       4. End the if statement and for loop
    9. Go to folder 'Plots' using cd('Plots') and make a folder in the name of current species the go into the spcecies folder.
    10. In the species folder open figure 1 and plot T vs Cp plot with linewidth as 3 and color red plot(T,Cp,'linewidth',3,'color','r') add proper label and titles to the plot and save it using saveas command as jpg file .
    11. Similarly open figure 2 and 3 and plot T vs H and T vs S with green and blue color respectively then after adding titles and label save them in current folder.
    12. Go back to the main grandparent directory in witch our program file and foler plots is saved .
    13. Make function molecular_weight for calculation of the weigth of the species .
        1. Take species name as input to the function .
        2. Define a matrix 'A' witch contain all the atoms of species written in file .
        3. Define a matrix 'W' witch contain weigth of atoms in A in respective position.
        4. Define a variable mw to zero .
        5. Create a for loop for each elements of matrix species
           1. Create a for loop for each elements of matrix 'A'
              1. Comapre each elements of A with elements of species using strcpmi built in funtion .
              2. If comparition is true add the weight from matrix 'W' to mw and save indice of A to a variable 'k'
           2. Check weather the elememt species is number by converting it to class double and camparing with number 1 beacuse str2num convert alphabets to an empty matrix
              
              1. If the if statements is varified add one minus the number times the weigth from the matrix 'W' in the positon k
    14. Print the molecular weigth of the species in command window with fprintf as fprintf('molecular weight of %s is %d\n',species,mw) and close the for loop
11. Close the flle with fclose

Main Program

close all
clear all
clc

% open THERMO.dat file for reading
f1=fopen('THERMO.dat','r');

%check for file opening error
if f1<0
    error('error opening file %s\n','THERMO.dat')
end

% Make a folder to save plots
mkdir('Plots')

%read line 1
line=fgetl(f1);

%read line 2 of global temp
global_temp=fgetl(f1);

%split line
A=strsplit(global_temp,' ');

%define global_temp
global_min_temp=str2num(A{2});
global_mid_temp=str2num(A{3});
global_max_temp=str2num(A{4});

%skipping 3  lines
fgetl(f1);
fgetl(f1);
fgetl(f1);


% for loop for 53 species
for i=1:53
    
    % LINE 1 of ith species
    line1=fgetl(f1);
    
    %split line 1
    L1=strsplit(line1,' ');
    
    %get species name and local temperature from line 1
    species = L1{1};
    local_min_temp = str2num(L1{end-3});
    local_max_temp = str2num(L1{end-2});
    local_mid_temp = str2num(L1{end-1});
    
    
    %LINE 2
    line2=fgetl(f1);
    
    %find position of E
    L2 = findstr(line2,'E');
    
    %get coefficient with the help of position of E
    a1 = str2num(line2(1:L2(1)+3));
    a2 = str2num(line2(L2(1)+4:L2(2)+3));
    a3 = str2num(line2(L2(2)+4:L2(3)+3));
    a4 = str2num(line2(L2(3)+4:L2(4)+3));
    a5 = str2num(line2(L2(4)+4:L2(5)+3));
    
    %line3
    line3=fgetl(f1);
    
    %find position of E
    L3 = findstr(line3,'E');
    
    %get coefficient with the help of position of E
    a6 = str2num(line3(1:L3(1)+3));
    a7 = str2num(line3(L3(1)+4:L3(2)+3));
    a8 = str2num(line3(L3(2)+4:L3(3)+3));
    a9 = str2num(line3(L3(3)+4:L3(4)+3));
    a10 = str2num(line3(L3(4)+4:L3(5)+3));
    
    %LINE 4
    line4=fgetl(f1);
    
    %find position of E
    L4 = findstr(line4,'E');
    
    %get coefficient with the help of position of E
    a11 = str2num(line4(1:L4(1)+3));
    a12 = str2num(line4(L4(1)+4:L4(2)+3));
    a13 = str2num(line4(L4(2)+4:L4(3)+3));
    a14 = str2num(line4(L4(3)+4:L4(4)+3));
    
    %Define temperature range
    T=linspace(local_min_temp,local_max_temp,500);
    
    % calculate Cp,H and S with the help of function
    [Cp,H,S] = thermodynamic_properties(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,...
        a11,a12,a13,a14,local_mid_temp,T);
    
    %go to folder plots
    cd('Plots')
    
    % make a folder for ith spceies
    mkdir(species)
    
    %go to species
    cd(species)
    
    %make plot of Cp
    figure(1)
    plot(T,Cp,'linewidth',3,'color','r');
    
    %label and title
    xlabel('Temperature(T)');
    ylabel('Specific heat(Cp)');
    title(['T vs Cp of ',species]);
    
    %save Cp plot
    saveas(figure(1),'T vs Cp.jpg');
    
    %make plot of H
    figure(2)
    plot(T,H,'linewidth',3,'color','g');
    
    %label and title
    xlabel('Temperature(T)');
    ylabel('Enthalpy(H)');
    title(['T vs H of ',species]);
    
    %save H plot
    saveas(figure(2),'T vs H.jpg');
    
    %make plot of S
    figure(3)
    plot(T,S,'linewidth',3,'color','b');
    
    %label and title
    xlabel('Temperature(T)');
    ylabel('Entropy(s)');
    title(['T vs S of ',species]);
    
    %save S plot
    saveas(figure(3),'T vs S.jpg')
    
    % go back to main directory
    cd('../..')
    
    %Calculate molecular weigth of the species
    mw = molecular_weight(species);
    
    %print molecular weigth of the species in command window
    fprintf('molecular weight of %s is %d\n',species,mw);
end

fclose(f1);

Functions :-

1) thermodynamic_properties

function [Cp,H,S ] = thermodynamic_properties(a1,a2,a3,a4,a5,a6,a7,a8,...
    a9,a10,a11,a12,a13,a14,local_mid_temp,T)


%Universal gas constant R
R = 8.314;

%For loop for calculating Cp,H and S for every value of temp range
for i=1:length(T)
    
    % HIGH-temperature coefficients
    if T(i)>=local_mid_temp
        Cp(i) = (a1 + a2*T(i) + a3*T(i)^2 + a4*T(i)^3 + a5*T(i)^4)*R;
        H(i) = (a1 + a2*T(i)/2 + a3*T(i)^2/3 + a4*T(i)^3/4 + a5*T(i)^4/5 + a6/T(i))*R*T(i);
        S(i) = (a1*log(T(i)) + a2*T(i) + a3*T(i)^2/2 + a4*T(i)^3/3 + a5*T(i)^4/4 + a7)*R;
        
    %LOW-temperature coefficients
    elseif T(i)<local_mid_temp
        Cp(i) = (a8 + a9*T(i) + a10*T(i)^2 + a11*T(i)^3 + a12*T(i)^4)*R;
        H(i) = (a8 + a9*T(i)/2 + a10*T(i)^2/3 + a11*T(i)^3/4 + a12*T(i)^4/5 + a13/T(i))*R*T(i);
        S(i) = (a8*log(T(i)) + a9*T(i) + a10*T(i)^2/2 + a11*T(i)^3/3 + a12*T(i)^4/4 + a14)*R;
    end
end

2) molecular_weight

function mw = molecular_weight(species)
%define a matrix for atoms
A = ['O', 'H', 'C', 'S', 'N' ,'A' ];

%define a matrix for corresponding weigths
W = [ 16 ,1 , 12 , 32, 14, 40];

% intial value to zero
mw = 0;

% for loop for speceis name
for i = 1:length(species)
    
    %nested for loop for camparing with A
    for j = 1:length(A)
        
        %Compare species with A without case sensitive
        if strcmpi(species(i),A(j))
            
            %add weigth to mw
            mw = mw + W(j);
            % save value to j to k
            k=j;
        end
    end
    
    %check for a number
    if str2num(species(i))>1
        
        % add weight of numtiple atoms to mw
        mw = mw + W(k)*(str2num(species(i))-1);
    end
end

RESULT

Screenshot of the window where all the folders are saved.

Command window

Plots for O2

1)

2)

3)

Plots for N2

1)

2)

3)

Plots for CO2

1)

2)

3)