Project 1 - Parsing NASA thermodynamic data

NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S using polynomials. They have also documented the co-efficients that are required to evaluate these polynomials for 1000+ species. The coefficients from the NASA thermodynamic file is extracted and the secific heat, enthalpy and entropy are calculated for different species and the results are plotted. I performed the following steps to write the code.

• First the given file format is converted as 'Therm.dat'.
• I used fopen command to open the given file format in the starting of the code .
• Using 'r' gave permission to read the file.
• After opening the file used fgetlcommand to read the first line of the file which is header, and the header is THERMO.
• Again fgetl command is used to read the second line of the file as trmperature range which contains global low, mid and high working temperatures.
• Next using strsplit command with delimeter ' ' (a space between inverted commas),command to to store the temperatures in cell array.
• str2num command is used to convert the temperatures from string to number.
• Temperature range is defined using linspace command. Which in this code takes 200 values between global low and high temperatures. Universal gas constant is defined R=8.314.
• The next 3 comments are skipped usingfgetl command which runs under a forloop.
• There are 53 species in the given file and for each species or element there are four lines, in which first line gives information about element name and number of atoms and temperature range. Second,third & fourth line contains 14 coefficients in which first seven are high temperature coefficients & the next seven are low temperature coefficients.fgetl command is used to read the 4 lines.
• strtok command is used to read the element name from the first line which takes the first word from that line.
• Next the character E is located using findstr command, which shows the position of the coefficients in that line. They are then stored in the variables a,b & c.
• As explained in the help video the coefficients are extracted. The low temperature coefficients are stored in the variables a1,a2,a3,a4,a5,a6,a7 and the high temperature coefficients are stored in the variables b1,b2,b3,b4,b5,b6,b7. the coefficients run for 53 times under a forloop.
• The coefficients for each elements get stored in the variables for each run.
• Separate functions for calculating specific heat, enthalpy, entropy, molecular weight and for plotting and they are called in the function for required action.
• The formulae for calculation is as follows.

The code is provided below.

clear all
close all
clc

%opening file and creating the header file
f=fopen('therm.dat','r');
Heading=fgetl(f);
temp_range=fgetl(f);

%assigning values to global low(glt),mid(gmt) and high(ght) temperatures 
gt=strsplit(temp_range,' ');
glt=str2num(gt{2});
gmt=str2num(gt{3});
ght=str2num(gt{4});

%defining working temperature range and universal gas constant R
T=linspace(glt,ght,200);
R=8.314;

%skipping the 3 comments
for i=1:3
    skip=fgetl(f)
end

for i=1:53
    %determining elements 
    l=fgetl(f);
    element=strtok(l);
    
    %reading next 3 lines which has 7 high and low temperature coefficients
    l1=fgetl(f);
    l2=fgetl(f);
    l3=fgetl(f);
    
    %locating E to extract the coefficients
    a=findstr(l1,'E');
    b=findstr(l2,'E');
    c=findstr(l3,'E');
    
    %storing the first 7 coefficients as high temperature coefficients
    a1=str2num(l1(1:a(1)+3));
    a2=str2num(l1(a(1)+4:a(2)+3));
    a3=str2num(l1(a(2)+4:a(3)+3));
    a4=str2num(l1(a(3)+4:a(4)+3));
    a5=str2num(l1(a(4)+4:a(5)+3));
    a6=str2num(l2(1:b(1)+3));
    a7=str2num(l2(b(1)+4:b(2)+3));
    
    %storing the first 7 coefficients as low temperature coefficients
    b1=str2num(l2(b(2)+4:b(3)+3));
    b2=str2num(l2(b(3)+4:b(4)+3));
    b3=str2num(l2(b(4)+4:b(5)+3));
    b4=str2num(l3(1:c(1)+3));
    b5=str2num(l3(c(1)+4:c(2)+3));
    b6=str2num(l3(c(2)+4:c(3)+3));
    b7=str2num(l3(c(3)+4:b(4)+3));
    
    %calculating enthalpy
    h=enthalpy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,T,R,glt,gmt,ght);
    
    %calculating entropy
    s=entropy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,T,R,glt,gmt,ght);
    
    %calculating specific heat
    cp=spec_heat(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,T,R,glt,gmt,ght);
    
    %plotting 
    plotting(cp,h,s,T,element);
    
    %calculating molecular weight of each species
    molecular_weight(i)=mol_weight(element);
end

1.Function for calculating the Specific heat.

For specific heat the input arguments required are low temperature coefficients, high temperature coefficients, global low temperature, global mid temperature, global high temperature, temperature range T and universal gas constant R.

Here in this code if else command is used to calculate the specific heat. If temperature range is less than or equal to global mid temperature low temperature coefficients are used to calculate the specific heat else high temperature coefficients are used.

function [cp]=spec_heat(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,T,R,glt,gmt,ght)

if T<gmt
    cp=R*(a1+a2.*T+a3.*(T.^2)+a4.*(T.^3)+a5.*(T.^4));
else     
    cp=R*(b1+b2.*T+b3.*(T.^2)+b4.*(T.^3)+b5.*(T.^4));
end

2.Function for calculating the enthalpy.

It requires same input arguments as that of specific heat.If temperature range is less than or equal to global mid temperature low temperature coefficients are used to calculate the enthalpy else high temperature coefficients are used.

function [h]=enthalpy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,T,R,glt,gmt,ght)

if T<gmt
    h=R*(a1+a2.*(T/2)+a3.*(T.^2)./3+a4.*(T.^3)./4+a5.*(T.^4)./5+a6./T);
else     
    h=R*(b1+b2.*(T/2)+b3.*(T.^2)./3+b4.*(T.^3)./4+b5.*(T.^4)./5+b6./T);
end

3.Function for calculating the entropy.

It requires same input arguments as that of specific heat.If temperature range is less than or equal to global mid temperature low temperature coefficients are used to calculate the entropy else high temperature coefficients are used.

function [s]=entropy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,T,R,glt,gmt,ght)
if T<=gmt
    s=R*(a1.*log(T)+a2.*T+a3.*(T.^2)./2+a4.*(T.^3)./3+a5.*(T.^4)./4+a7);
else
    s=R*(b1.*log(T)+b2.*T+b3.*(T.^2)./2+b4.*(T.^3)./3+b5.*(T.^4)./4+b7);
end

4.Plotting

For plotting there is no any output argument and the input arguments are specific heat, enthalpy, entropy, temperature range T and the element name.

• Used mkdir command to create a new folder in the current current folder (pwd) and saved it in the name of individual elements.
• The folders of each element with its name will be created as the code runs for 53 times.
• After that used appropriate variables for plotting the graph with the specific element name,with appropriate x & y labels, and then used saveas command to save all the three graphs with specific names and file format in the folders created.
• Then used " cd .. " to change the folder to one folder up.

function plotting(cp,h,s,T,element)
cd(pwd)
dir_name=element;
mkdir(dir_name);

figure(1)
plot(T,cp,'lineWidth',3)
title(sprintf('Specific Heat vs Temperature rn for %s:',element))
xlabel('Temperature')
ylabel('Specific Heat')

figure(2)
plot(T,h,'linewidth',3)
title(sprintf('Enthalpy vs Temperature rn for %s:',element))
xlabel('Temperature')
ylabel('Enthalpy')

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

cd(dir_name)
saveas(1,'cp_vs_T.png')
saveas(2,'h_vs_T.png')
saveas(3,'s_vs_T.png')

cd ..

end

5.Molecular weight

This function for molecular weight calculations requires element name as input argument and it is saved in the variable 'element' and the output argument is molecular weight.

• 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 repeatition of any element, str2num is used to find that and (n-1) is multiplied to atomic_wt of element to calculate molecular weight.

function [w] = mol_weight(element)


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

The plots generated for the following element is as follows

1.O2

2.N2

3.CO2

References used for the projects is

1.Google.

2.Skill-lync project website.

3.Help forum of matlab.