PARSING NASA THERMODYNAMIC DATA AND CALCULATING THERMAL PROPERTIES OF VARIOUS GAS SPECIES FROM THE DATA USING MATLAB

AIM :

To parse the NASA thermodynamic data from the data file and to calculate the thermodynamic properties like Specific heat at constant pressure(Cp), Enthalpy (H)and Entropy(S) of various species present in the data file using MATLAB programming.

Data File :

https://drive.google.com/file/d/1AikCHDlz9s_fu8whE81qee_GD4AyYpCP/view?usp=sharing

NASA came up up with the polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S. NASA documented the coefficients that are required to evaluate these polynomials for 1000+ species.

This a sample data for a single species .The 2nd line indicates the Global temperature ranges and the 6th line indicates the species name and its local temperature ranges and the next three lines indicate the co-efficients that have to be evaluated for finding the thermodynamic properties.

For each species the data file contains 14 coefficents

Cp/R = a1 + a2 T + a3 T² + a4 T³ + a5 T⁴

H/RT = a1+ a2 T/2 + a3 T² /3+ a4 T³ /4+ a5 T⁴/5 + a6

S/R = a1 lnT + a2 T + a3 T² /2 + a4 T³ /3 + a5 T⁴/4 + a7

Cp/R = a8 + a9 T + a10 T² + a11 T³ + a12 T⁴

H/RT = a8 + a9 T/2 + a10 T² /3+ a11 T³ /4+ a12 T⁴/5 + a13

S/R = a8 lnT + a9 T + a10 T² /2 + a11 T³ /3 + a12 T⁴/4 + a14

a1 to a7 are the low temperature coefficients and a8 to a14 are high temperature coefficients

a1 to a14 values are calculated from the NASA data file

File parsing :

It is the method of splitting the file or other input in to a piece of data that can be easily stored or manipulated. The NASA data contains 53 gas species and total 213 lines of coefficient data.

MATLAB CODE :

clear all
close all
clc
%opening and reading the NASA Thermo data file
file=fopen('THERMO.dat','r')
%heading
line1=fgetl(file);
%Global temperature values line
line2=fgetl(file);
%splitting the total string in to cell arrays
A=strsplit(line2,' ');
glowtemp=str2double(A{2}); gmidtemp=str2double(A{3}); ghightemp=str2double(A{4});
%reading the comments in the line
for z=1:3
    fgetl(file);
end
%Getting in to the species block and total species are 53
a=strings([1,16]); c=[]; %taking empty string array and empty vector array for storing the coefficients
R_uni=8.314;%Universal gas constant in (Kj/Kmol-K)
species_w=fopen('species_mass.txt','w')%creating a text file to store molecular weights of species
for i=1:53
    speciesline=fgetl(file);%species block title line
    %Extracting the species name
    speciesname=speciesline(1:10); %By using findstr and size commands find the position of charecters
    speciesname=deblank(speciesname);%Removes blank spaces 
    %Extracting Local temperatures
    speciestemp=speciesline(49:75);%This is the location where temperature values are present.The positions can be known by findstr command
    localtempstring=strsplit(speciestemp,' ');%splitting the local temperature string into cell array
    locallowtemp=str2double(localtempstring{1});
    localmidtemp=str2double(localtempstring{3});
    localhightemp=str2double(localtempstring{2});
    localtempvalues=[locallowtemp localmidtemp localhightemp];%storing local temp values in an array
    temprange1=linspace(locallowtemp,localmidtemp,500);
    temprange2=linspace(localmidtemp,localhightemp,500);%range of local temp values
    coeff=fgetl(file);%species block coefficients 1st line
    coeff1=findstr(coeff,'E');
    %Extracting the coefficients
    a(1)=coeff(1:coeff1(1)+3);
    c(1)=str2double(a(1));
    for m=1:4
        a(m+1)=coeff(coeff1(m)+4:coeff1(m+1)+3);
        c(m+1)=str2double(a(m+1));
    end
    newline=fgetl(file);%species block coefficients 2nd line
    coeff2=findstr(newline,'E');
    a(6)=newline(1:coeff2(1)+3);
    c(6)=str2double(a(6));
    for p=1:4
        a(p+6)=newline(coeff2(p)+4:coeff2(p+1)+3);
        c(p+6)=str2double(a(p+6));
    end
    nextnewline=fgetl(file);%species block coefficients 3rd line
    coeff3=findstr(nextnewline,'E');
    a(11)=nextnewline(1:coeff3(1)+3);
    c(11)=str2double(a(11));
    for q=1:3
         a(q+11)=nextnewline(coeff3(q)+4:coeff3(q+1)+3);
        c(q+11)=str2double(a(q+11));
    end
a1=c(1); a2=c(2); a3=c(3); a4=c(4); a5=c(5); a6=c(6); a7=c(7); a8=c(8); a9=c(9); a10=c(10); a11=c(11); a12=c(12); a13=c(13); a14=c(14);
%calculation of species molecular weight
molweight=weight(speciesname);
%charecterstic gas constant in (Kj/Kg-K)
R=(R_uni/molweight)
%Evaluation of specific heat , Enthalpy and Entropy
[cp1 H1 S1 cp2 H2 S2]=evaluateval(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,temprange1,temprange2,R);
cur_dir=pwd;
cd('C:\Users\hp\Documents\MATLAB')%locating the directory to create species folders
mkdir(speciesname);
cd(speciesname);
%plotting and saving
figure(1)
plot(temprange1,cp1,'linewidth',3,'color','r')
hold on
plot(temprange2,cp2,'linewidth',3,'color','g')
hold off
xlabel('temperature (K)')
ylabel('specific heat (Kj/kg-k)')
title(sprintf('specific heat vs temperature range of %s',speciesname))
legend('low temp','High temp','location','southeast')
xlim([0 5000])
xticks(0:500:5000)
saveas(figure(1),'specific heat.jpg')
figure(2)
plot(temprange1,H1,'linewidth',3,'color','b')
hold on
plot(temprange2,H2,'linewidth',3,'color','y')
hold off
xlabel('temperature (K)')
ylabel('Enthalpy (Kj/kg)')
title(sprintf('Enthalpy vs temperature range of %s',speciesname))
legend('low temp','High temp','location','southeast')
xlim([0 5000])
xticks(0:500:5000)
saveas(figure(2),'Enthalpy.jpg')
figure(3)
plot(temprange1,S1,'linewidth',3,'color','m')
hold on
plot(temprange2,S2,'linewidth',3,'color','c')
hold off
xlabel('temperature (K)')
ylabel('Entropy (Kj/kg)')
title(sprintf('Entropy vs temperature range of %s',speciesname))
legend('low temp','High temp','location','southeast')
xlim([0 5000])
xticks(0:500:5000)
saveas(figure(3),'Entropy.jpg')
cd(cur_dir);
%storing molecular weights in text file
fprintf(species_w,'molecular mass of %s is %f',speciesname,molweight)
fprintf(species_w,'\n')
end
fclose(species_w)
fclose(file)

FUNCTION CODE for evaluating Cp, H and S :

function [cp1 H1 S1 cp2 H2 S2] = evaluateval(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,temprange1,temprange2,R)
for j=1:length(temprange1)
        cp1(j)=R*(a8+(a9*temprange1(j))+(a10*(temprange1(j)).^2)+(a11*(temprange1(j)).^3)+(a12*(temprange1(j)).^4));
         H1(j)=R*temprange1(j)*(a8+(a9*(temprange1(j)/2))+(a10*((temprange1(j)).^2)/3)+(a11*((temprange1(j)).^3/4))+(a12*((temprange1(j)).^4/5))+a13/temprange1(j));
         S1(j)=R*((a8*log(temprange1(j)))+(a9*(temprange1(j)))+(a10*((temprange1(j)).^2)/2)+(a11*((temprange1(j)).^3/3))+(a12*((temprange1(j)).^4/4))+a14);
end
for w=1:length(temprange2)
    cp2(w)=R*(a1+(a2*temprange2(w))+(a3*(temprange2(w)).^2)+(a4*(temprange2(w)).^3)+(a5*(temprange2(w)).^4));
         H2(w)=R*temprange2(w)*(a1+(a2*(temprange2(w)/2))+(a3*((temprange2(w)).^2)/3)+(a4*((temprange2(w)).^3/4))+(a5*((temprange2(w)).^4/5))+a6/temprange2(w));
         S2(w)=R*((a1*log(temprange2(w)))+(a2*(temprange2(w)))+(a3*((temprange2(w)).^2)/2)+(a4*((temprange2(w)).^3/3))+(a5*((temprange2(w)).^4/4))+a7);
end

FUNCTION CODE for calculation of Species Molecular weight :

function molweight= weight(speciesname)
element={'H','O','N','C','S'};
atomicmass=[1.008 15.999 14.007 12.011 32.065];
molweight=0;
if (strcmp(speciesname,'AR'))==0
for i=1:length(speciesname)
    for j=1:length(element)
        if strcmp(speciesname(i),element(j))
            molweight=molweight+atomicmass(j);
            x=j;
          end
    end
    y=str2num(speciesname(i));
    if y>1
        molweight=molweight+atomicmass(x)*(y-1);
    end
end
else
    molweight=39.908;
end

Explanation of Code in stepwise manner :

1) First the file has to be read for that open the THERMO.dat file using file open command

2) fgetl command reads the first line of the data file if that command is again executed then it reads the 2nd line . So this command has to be written every time for reading each and every line in data file.

3)The 2nd line of data file contains the global temperature values. Using the strsplit command the entire string is broken into cell arrays and the cell string arrays are converted into numerical values using str2double command

4) The next comment lines can be skipped by using the fgetl command .

5) The Header data is succesfully read and now the species data has to be read. For reading the entire 53 species data a For loop is created.

6) Create an Empty string array and empty data array for storing the coefficients in the form of strings and values.

7) By using findstr and size commands find the position of charecters in the first line of species and now the species name and the species local temperatures ranges are to be extracted.

8)For reading the species name first 18 charecters can be taken and the local temperatures are present from the charecter position 49 to 75

9) Store the local temperature ranges from low to high and create a array for storing the local temperature values and also define a low temperature range and a high temperature range for evaluating Cp1,H1,S1 and Cp2, H2 and S2.

10) All the 14 coefficients can be calculated bu using different nested for loops

11) The Specific heat, Enthalpy and Entropy for lower and higher temperatures are calculated by taking the seperate function called evaluate and calling the function into the main function.

12) Similarly molecular weight of each species has to be calculated in a seperate fuction called weight for calculating the charecterstic gas constant(R) of individual species (R = R_universal/molecular weight)

13) plotting all the values of Cp, H and S wrt to high and low temperatures.

14) Creating a separate text file(species_mass.txt) for writing the molecular weight values of all the species.

15)creating seperate folder for each species and saving the plots in it

16)Terminating the final loop and closing all the files.

Outputs :

Similarly for every 53 species the above three plots are generated.

All SPECIES OUTPUT PLOTS IN FOLDERS :

https://drive.google.com/drive/folders/1q28J5cbBCFy4S0o-4j3vmFFycOs5P01S?usp=sharing

Molecular weight of all species in text file :

https://drive.google.com/file/d/1aOOE5hLzisultrzDg3gyRFtWBMP9GnXY/view?usp=sharing

Errors:

All are typo errors and they bare rectified using the error message in command propmt but proper concentration is required for writing the code for molecular weight function the code doesnt worked for calculating AR molecular weight and it can be rectified by taking the special if condition.