Challenge

Project 1 - Parsing NASA thermodynamic data

KKiran More9 Mar 2022

Aim:

To write MATLAB CODE to extract data from NASA thermodynamic data file.
To calculate thermodynamic properties like specific heat, enthalpy and entropy for all species present in the data file, from data given in file and to plot variation of specific heat, enthalpy, and entropy of each species with respect to temperature.
To calculate molecular weight of each species given in a data file.

Introduction:

File parsing means reading data from given file to facilitate some transformation on the data present in the file. In the current project parsing of NASA thermodynamic data is done using MATLAB code. File is present in .dat format. File contains 53 species. Each species has 14 coefficients by using which specific heat, enthalpy and entropy are calculated. Out of 14 coefficients, first 7 are high-temperature coefficients and next 7 are low-temperature coefficients. Also, molecular weight of each species is calculated.

After calculating specific heat, enthalpy and entropy, variation of it is plotted with respect to temperature for each species. Each species will have three plots. So total number of plots are 53x3=159. In the code arrangement is done in a such way that, all plots are saved in single folder species wise separately.

Following mathematical relations are used to calculate specific heat, enthalpy, and entropy, respectively.

Where Cp is specific heat at constant pressure, R is universal gas constant, H is enthalpy, S is entropy and T is temperature.

Code explanation

Given .dat file is read by using ‘fopen(' THERMO.dat', 'r')
```
f1 = fopen('THERMO.dat', 'r');
```
By using ‘fget(l)’ command, first line (THERMO)is read.
```
fgetl(f1);
```
Global temperatures (2^nd line of file) is read by using ‘aa = fgetl(f1)’. Then this line is spitted by using ‘fget(l)’ command. Then global temperatures (300, 1000 and 5000) are read.
```
aa = fgetl(f1);
tt = strsplit(aa, ' ');
global_low_temp = str2num(tt{2});
global_mid_temp = str2num(tt{3});
global_high_temp = str2num(tt{4});
```

Line no 3, 4 and 5 are read using ‘fget(l)’ command.

 To eliminate comment lines
for i=3:5
     fgetl(f1);
end

Now each species has 4 lines. First line contains the species name and local temperatures (low, middle and high temperature). Next line
(O L 1/90O 1 G 200.000 3500.000 1000.000 1)
is read by using ‘fget(l)’ command. Then to identify the species, above line is splitted using ‘fget(l)’ command and name species is identified.
```
for j=1:53
tline=fgetl(f1);
A = strsplit(tline, ' ');
species = (A{1});
```

Now to determine the values of low, middle, and high temperature, the location of dot(.) is identified and with respect to position of dot(.) values of low, middle, and high temperature are read.

pos_dot=strfind(tline,'.');
low_temp = tline(pos_dot(1)-3:pos_dot(1)+3);
low_temp = str2num(low_temp);
high_temp = tline(pos_dot(2)-4:pos_dot(2)+3);
high_temp =  str2num(high_temp);
mid_temp = tline(pos_dot(3)-4:pos_dot(3)+3);
mid_temp =  str2num(mid_temp);

Similarly, next 3 lines are read using ‘fget(l)’ command. First and second line has 5 coefficients each and 3^rd line has 4 coefficients. To read value of each coefficient, position of ‘E’ has been found out in each line and with respect to position of ‘E’ each coefficient value is read.
```
dataline_1 = fgetl(f1);
dataline_2 = fgetl(f1);
dataline_3 = fgetl(f1);

A = strfind(dataline_1,'E');
B = strfind(dataline_2,'E');
C = strfind(dataline_3,'E');
```
Temperature values are specified using ‘linspace’ command and value of R is defined.
```
T = linspace(low_temp, high_temp,500);
R = 8.314;
```
Specific heat, enthalpy and Entropy is calculated.

for i=1:length(T)
    if T(i)>mid_temp
        cp(i)= (t1+t2.*T(i)+t3.*T(i).^2+t4.*T(i).^3+t5.*T(i).^4)*R;
        H(i)= R*T(i).*(t1+((t2*T(i)/2))+((t3*T(i).^2)/3)+((t4*T(i).^3)/4)+((t5*T(i).^4)/5)+((t6./T(i))));
        S(i)= (t1.*log(T(i))+t2.*T(i)+t3.*(T(i)^2)/2+t4.*(T(i)^3)/3+t5.*(T(i)^4)/4+t7)*R;
    else
        cp(i)= (t8+t9.*T(i)+t10.*T(i).^2+t11.*T(i).^3+t12.*T(i).^4)*R;
        H(i)= R*T(i).*(t8+((t9*T(i)/2))+((t10*T(i).^2)/3)+((t11*T(i).^3)/4)+((t12*T(i).^4)/5)+((t13./T(i))));
        S(i)= (t8.*log(T(i))+t9.*T(i)+t10.*(T(i)^2)/2+t11.*(T(i)^3)/3+t12.*(T(i)^4)/4+t14)*R;
    end
end

10. Then variation of specific heat, enthalpy and entropy are plotted with respect to temperature.

% To plot varation of cp with T
figure(1)
plot(T, cp, 'linewidth',2,'color','g')
xlabel('Temperature in K')
ylabel('Specific heat in KJ/kg.K')
grid on
title(sprintf('Variation of specific heat with temperature for %s', species))
saveas(figure(1),'cp.png');

% To plot variation of enthapy with temperarure
figure(2)
plot(T, H, 'linewidth',2,'color','r')
xlabel('Temperature in K')
ylabel('Enthalpy in KJ/kg')
grid on
title(sprintf('Variation of Enthalpy with temperature for %s', species))
saveas(figure(2),'H.png');

% To plot variation of entropy with temperarure
figure(3)
plot(T, S, 'linewidth',2,'color','b')
xlabel('Temperature in K')
ylabel('Entropy in KJ/K')
grid on
title(sprintf('Variation of Entropy with temperature for %s', species))
saveas(figure(3),'S.png');

11. Molecular weight is calculated.

% To calculate molecular weight
 Molecular_wt = molecular_weight(species);
 k = sprintf('Molecular weight of %s is %f', species,Molecular_wt);
 disp(k)

12. One folder is created named ‘species’ and all plots are saved in that folder, species wise separately.

% Saving in a folder
cd('C:\Users\Yogesh\OneDrive\Documents\MATLAB\species')
mkdir(species);
cd(species);

Function code to calculate Molecular weight

function[molecular_weight] = molecular_weight(species)
elements = {'H','C','O','N','A','S'};
atomic_weight = [1,12,16,14,39,32];
molecular_weight=0;
 for i=1:length(species)
     for j=1:length(elements)
         if strcmp(species(i),elements(j))
             molecular_weight = molecular_weight+atomic_weight(j) 
             a=j;
         end
         end  
       n = str2num(species(i));
       if(n>1)
             molecular_weight = molecular_weight+(atomic_weight(a)*(n-1))
     end
 end
end

Complete Code

clear all
close all
clc

% Getting data from provided file
f1 = fopen('THERMO.dat', 'r');
fgetl(f1);
aa = fgetl(f1);
tt = strsplit(aa, ' ');
global_low_temp = str2num(tt{2});
global_mid_temp = str2num(tt{3});
global_high_temp = str2num(tt{4});

% To eliminate comment lines
for i=3:5
     fgets(f1)
end

% To get species data 
for j=1:53
    tline=fgetl(f1);
    A = strsplit(tline, ' ');
    species = (A{1});
    pos_dot=strfind(tline,'.');
    low_temp = tline(pos_dot(1)-3:pos_dot(1)+3);
    low_temp = str2num(low_temp);
    high_temp = tline(pos_dot(2)-4:pos_dot(2)+3);
    high_temp =  str2num(high_temp);
    mid_temp = tline(pos_dot(3)-4:pos_dot(3)+3);
    mid_temp =  str2num(mid_temp);
    
    dataline_1 = fgetl(f1);
    dataline_2 = fgetl(f1);
    dataline_3 = fgetl(f1);

    A = strfind(dataline_1,'E');
    B = strfind(dataline_2,'E');
    C = strfind(dataline_3,'E');

    t1 = dataline_1(1:A(1)+3);
    t1 = str2num(t1);
    t2 = dataline_1(A(1)+4:A(2)+3);
    t2 = str2num(t2);
    t3 = dataline_1(A(2)+4:A(3)+3);
    t3 = str2num(t3);
    t4 = dataline_1(A(3)+4:A(4)+3);
    t4 = str2num(t4);
    t5 = dataline_1(A(4)+4:A(5)+3);
    t5 = str2num(t5);

    t6 = dataline_2(1:B(1)+3);
    t6 = str2num(t6);
    t7 = dataline_2(B(1)+4:B(2)+3);
    t7 = str2num(t7);
    t8 = dataline_2(B(2)+4:B(3)+3);
    t8 = str2num(t8);
    t9 = dataline_2(B(3)+4:B(4)+3);
    t9 = str2num(t9);
   t10 = dataline_2(B(4)+4:B(5)+3);
   t10 = str2num(t10);

   t11 = dataline_3(1:C(1)+3);
   t11 = str2num(t11);
   t12 = dataline_3(C(1)+4:C(2)+3);
   t12 = str2num(t12);
   t13 = dataline_3(C(2)+4:C(3)+3);
   t13 = str2num(t13);
   t14 = dataline_3(C(3)+4:C(4)+3);
   t14 = str2num(t14);


T = linspace(low_temp, high_temp,500);
R = 8.314;

% specific heat, Entahlpy and Entropy calculation
for i=1:length(T)
    if T(i)>mid_temp
        cp(i)= (t1+t2.*T(i)+t3.*T(i).^2+t4.*T(i).^3+t5.*T(i).^4)*R;
        H(i)= R*T(i).*(t1+((t2*T(i)/2))+((t3*T(i).^2)/3)+((t4*T(i).^3)/4)+((t5*T(i).^4)/5)+((t6./T(i))));
        S(i)= (t1.*log(T(i))+t2.*T(i)+t3.*(T(i)^2)/2+t4.*(T(i)^3)/3+t5.*(T(i)^4)/4+t7)*R;
    else
        cp(i)= (t8+t9.*T(i)+t10.*T(i).^2+t11.*T(i).^3+t12.*T(i).^4)*R;
        H(i)= R*T(i).*(t8+((t9*T(i)/2))+((t10*T(i).^2)/3)+((t11*T(i).^3)/4)+((t12*T(i).^4)/5)+((t13./T(i))));
        S(i)= (t8.*log(T(i))+t9.*T(i)+t10.*(T(i)^2)/2+t11.*(T(i)^3)/3+t12.*(T(i)^4)/4+t14)*R;
    end
end

% Saving in a folder
cd('C:\Users\DELL\Documents\MATLAB\Parsing')
mkdir(species);
cd(species);

% To plot varation of cp with T
figure(1)
plot(T, cp, 'linewidth',2,'color','g')
xlabel('Temperature in K')
ylabel('Specific heat in KJ/kg.K')
grid on
title(sprintf('Variation of specific heat with temperature for %s', species))
saveas(figure(1),'cp.png');

% To plot variation of enthapy with temperarure
figure(2)
plot(T, H, 'linewidth',2,'color','r')
xlabel('Temperature in K')
ylabel('Enthalpy in KJ/kg')
grid on
title(sprintf('Variation of Enthalpy with temperature for %s', species))
saveas(figure(2),'H.png');

% To plot variation of entropy with temperarure
figure(3)
plot(T, S, 'linewidth',2,'color','b')
xlabel('Temperature in K')
ylabel('Entropy in KJ/K')
grid on
title(sprintf('Variation of Entropy with temperature for %s', species))
saveas(figure(3),'S.png');
 
cd('C:\Users\DELL\Documents\MATLAB\Parsing')
 
% To calculate molecular weight
 Molecular_wt = molecular_weight(species);
 k = sprintf('Molecular weight of %s is %f', species,Molecular_wt);
 disp(k)
end

Results:

Coeficients

t1 =

    4.0172


t2 =

    0.0022


t3 =

  -6.3366e-07


t4 =

   1.1425e-10


t5 =

  -1.0791e-14


t6 =

  111.8567


t7 =

    3.7851


t8 =

    4.3018


t9 =

   -0.0047


t10 =

   2.1158e-05


t11 =

  -2.4276e-08


t12 =

   9.2923e-12


t13 =

  294.8080


t14 =

    3.7167