Parsing NASA Thermodynamic data using matlab

AIM:

         To write a code in Matlab to parse the NASA thermodynamic data file and then calculate thermodynamic properties of various gas species.

OBJECTIVE:

  • To write a function that extracts the 14 co-efficients of each species in the data file.
  • To calculate the enthalpy, entropy and specific heats for all the species in the data file.
  • To calculate the molecular weight of each species and display it in the command window.
  • To plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species and save them.

INTRODUCTION:

 

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

        This is done by parsing the NASA thermodynamic data file 'THERMO.dat' to extract the 14 coefficients of each species to calculate thermodynamic properties such as Cp,H and S by using the below formulas.

                1

         Here,

              a1,a2,a3,a4,a5,a6,a7 are coefficients

              T is temperature in kelvin

              R is universal gas constant

              Cp is specific heat (J/mol.K)

              H is enthalpy (J)

              S is entropy (J/K)

             2

A VIEW OF NASA THERMODYNAMIC FILE:

           3

Here,

  • The 1st five line represents the header
  • The line next to the header contains the individual species data such as local temperature and numerical coefficients for NASA polynomials which is unique for each species.
  •  Here, in the local temperature 

                                200 represents low local temperature.

                                3500 represents local high temperature.

                                1000 represents local mid temperature.

  • In numerical coefficients for NASA polynomials the FIRST 7 coefficients represent HIGH-temperature coefficients and the SECOND 7 coefficients represent LOW-temperature coefficients.

MAIN PROGRAM:

clear all
close all
clc

% Opening the thermo.dat file
f1 = fopen('thermo.dat','r');

% Skipping first 5 lines
for i=1:5
    skippinglines=fgetl(f1);
end

% Defining the value of R(J/mol.K)
R=8.3145;

% Function that extracts the 14 co-efficients
for i=1:53
    
% Extracting line 1 
line1=fgetl(f1);
[name,local_temp]=strtok(line1,'G');

% Extracting name of species
species_name=strtok(name);

% Extracting the local temperature
temp=strsplit(local_temp);
local_low_temp=str2double(temp{2});
local_high_temp=str2double(temp{3});
local_mid_temp=str2double(temp{4});
T=linspace(local_low_temp,local_high_temp,350);


% Extracting the higher and lower temperature coefficients
line2=fgetl(f1);
line3=fgetl(f1);
line4=fgetl(f1);

a=strfind(line2,'E');
b=strfind(line3,'E');
c=strfind(line4,'E');

% Higher temperature coefficients
a1=str2double(line2(1:a(1)+3));
a2=str2double(line2(a(1)+4:a(2)+3));
a3=str2double(line2(a(2)+4:a(3)+3));
a4=str2double(line2(a(3)+4:a(4)+3));
a5=str2double(line2(a(4)+4:a(5)+3));
a6=str2double(line3(1:b(1)+3));
a7=str2double(line3(b(1)+4:b(2)+3));

% Lower temperature coefficients
b1=str2double(line3(b(2)+4:b(3)+3));
b2=str2double(line3(b(3)+4:b(4)+3));
b3=str2double(line3(b(4)+4:b(5)+3));
b4=str2double(line4(1:c(1)+3));
b5=str2double(line4(c(1)+4:c(2)+3));
b6=str2double(line4(c(2)+4:c(3)+3));
b7=str2double(line4(c(3)+4:c(4)+3));

% Calculating specific heat for all species
Cp=specific_heat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,local_mid_temp);

% Calculating enthalpy
H=enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,R,T,local_mid_temp);

%Calculating entropy
S=entropy(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,R,T,local_mid_temp);

% Calculating molecular weight
atoms=['H','C','O','N','R','S'];
atomic_weight=[1 12 16 14 40 32];
MW=0;
for i=1:length(species_name)
    for j=1:length(atoms)
        if strcmp(species_name(i),atoms(j))
            MW=MW+atomic_weight(j);
            dt=j;
        end
    end
    n=str2double(species_name(i));
    if n>1
        MW=MW+(atomic_weight(dt)*(n-1));
    end
end
fprintf('Molecular weight of %s:%f',species_name,MW)
disp(' ');

% Creating file directory
mkdir(['E:matlab filesproject7',species_name]);
cd(['E:matlab filesproject7',species_name]);

%plotting specific heat vs temperature
figure(1)
plot(T,Cp,'linewidth',2,'color','b')
xlabel('TEMPERATURE')
ylabel('SPECIFIC HEAT')
saveas(gcf,'TEMPERATURE VS SPECIFIC HEAT.png')

%plotting enthalpy vs temperature
figure(2)
plot(T,H,'linewidth',2,'color','g')
xlabel('TEMPERATURE')
ylabel('ENTHALPY')
saveas(gcf,'TEMPERATURE VS ENTHALPY.png')


%plotting entropy vs temperature
figure(3)
plot(T,S,'linewidth',2,'color','r')
xlabel('TEMPERATURE')
ylabel('ENTROPY')
saveas(gcf,'TEMPERATURE VS ENTROPY.png')
cd(pwd)

file='E:matlab filesproject7';
cd(file)

end

PROGRAM FOR SPECIFIC HEAT FUNCTION:

% Specific heat function
function[Cp]=specific_heat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,local_mid_temp)
if T>local_mid_temp
    
    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

PROGRAM FOR ENTHALPY FUNCTION:

% Enthalpy function
function[H]=enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,R,T,local_mid_temp)
if T>local_mid_temp
    H=R.*T.*(a1+(a2.*T)./2+(a3.*T.^2)./3+(a4.*T.^3)./4+(a5.*T.^4)./5+a6./T);
else
    H=R.*T.*(b1+(b2.*T)./2+(b3.*T.^2)./3+(b4.*T.^3)./4+(b5.*T.^4)./5+b6./T);
end

PROGRAM FOR ENTROPY FUNCTION:

% Entropy function
function[S]=entropy(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,R,T,local_mid_temp)
if T>local_mid_temp
    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

EXPLANATION OF PROGRAM:

STEP1:

clear all
close all
clc

% Opening the thermo.dat file
f1 = fopen('thermo.dat','r');

% Skipping first 5 lines
for i=1:5
    skippinglines=fgetl(f1);
end

% Defining the value of R(J/mol.K)
R=8.3145;
  • In the above program the 1st three lines are used to clear the previous files data in both command window and in the workspace.
  • Then using the fopen command we are reading the NASA thermodynamic data file.
  • Then using the for command we are skipping the first 5 lines of data file at it contains comments and unwanted datas.
  • Since we know that R represents universal gas constant hence we have defined its value.

STEP2:

% Function that extracts the 14 co-efficients
for i=1:53
    
% Extracting line 1 
line1=fgetl(f1);
[name,local_temp]=strtok(line1,'G');

% Extracting name of species
species_name=strtok(name);
  • In the above code we have used for loop to represent all the species in the data file
  • Then we are extracting the line 1 using strtok command to seperate the species name and the temperature using a delimiter 'G'
  • Then we are using the strtok command over the name variable to extract the species name

STEP3:

% Extracting the local temperature
temp=strsplit(local_temp);
local_low_temp=str2double(temp{2});
local_high_temp=str2double(temp{3});
local_mid_temp=str2double(temp{4});
T=linspace(local_low_temp,local_high_temp,350);
  • In the above code we have used the strsplit command over the variable local_temp to split the variable such that we can access individual temperatures and provide each temperature a seperate variable name.
  • To define the temperature from low to high over a range of 350 values we have used the linspace command.

STEP4:

% Extracting the higher and lower temperature coefficients
line2=fgetl(f1);
line3=fgetl(f1);
line4=fgetl(f1);

a=strfind(line2,'E');
b=strfind(line3,'E');
c=strfind(line4,'E');

% Higher temperature coefficients
a1=str2double(line2(1:a(1)+3));
a2=str2double(line2(a(1)+4:a(2)+3));
a3=str2double(line2(a(2)+4:a(3)+3));
a4=str2double(line2(a(3)+4:a(4)+3));
a5=str2double(line2(a(4)+4:a(5)+3));
a6=str2double(line3(1:b(1)+3));
a7=str2double(line3(b(1)+4:b(2)+3));

% Lower temperature coefficients
b1=str2double(line3(b(2)+4:b(3)+3));
b2=str2double(line3(b(3)+4:b(4)+3));
b3=str2double(line3(b(4)+4:b(5)+3));
b4=str2double(line4(1:c(1)+3));
b5=str2double(line4(c(1)+4:c(2)+3));
b6=str2double(line4(c(2)+4:c(3)+3));
b7=str2double(line4(c(3)+4:c(4)+3));
  • In the above code we are using fgetl command to extract the coefficients in line 2,3, and 4.
  • Then we are using strfind command to locate letter 'E' in each line.
  • Then the strings are converted to values using str2double command
  • Then we are assigning a1,a2,a3,a4,a5,a6,a7 as variable name for first 7 coefficients and b1,b2,b3,b4,b5,b6,b7 for next 7 coefficients using the relation.

STEP5:

% Calculating specific heat for all species
Cp=specific_heat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,local_mid_temp);

% Calculating enthalpy
H=enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,R,T,local_mid_temp);

%Calculating entropy
S=entropy(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,R,T,local_mid_temp);

% Calculating molecular weight
atoms=['H','C','O','N','R','S'];
atomic_weight=[1 12 16 14 40 32];
MW=0;
for i=1:length(species_name)
    for j=1:length(atoms)
        if strcmp(species_name(i),atoms(j))
            MW=MW+atomic_weight(j);
            dt=j;
        end
    end
    n=str2double(species_name(i));
    if n>1
        MW=MW+(atomic_weight(dt)*(n-1));
    end
end
fprintf('Molecular weight of %s:%f',species_name,MW)
disp(' ');
  • In the above program we are calling in the function of specific heat, enthalpy and entropy to calculate the values for all the species in the data file
  • The molecular weight is also calculated for all the species using the for loop.

STEP6:

% Creating file directory
mkdir(['E:matlab filesproject7',species_name]);
cd(['E:matlab filesproject7',species_name]);

%plotting specific heat vs temperature
figure(1)
plot(T,Cp,'linewidth',2,'color','b')
xlabel('TEMPERATURE')
ylabel('SPECIFIC HEAT')
saveas(gcf,'TEMPERATURE VS SPECIFIC HEAT.png')

%plotting enthalpy vs temperature
figure(2)
plot(T,H,'linewidth',2,'color','g')
xlabel('TEMPERATURE')
ylabel('ENTHALPY')
saveas(gcf,'TEMPERATURE VS ENTHALPY.png')


%plotting entropy vs temperature
figure(3)
plot(T,S,'linewidth',2,'color','r')
xlabel('TEMPERATURE')
ylabel('ENTROPY')
saveas(gcf,'TEMPERATURE VS ENTROPY.png')
cd(pwd)

file='E:matlab filesproject7';
cd(file)

end

 

  • In the above code we have used mkdir command which creates the directory dirname in the current directory, if dirname represents a relative path to save the plots of each species seperately.
  • The cd command is used to display the new current folder.
  • Then finally we are ending the for loop using end command.

Screenshot of the window where all the folders are saved.

       5

PLOTS FOR o2:

T Vs H

                       3

T Vs S

                    

T Vs Cp

                   

PLOTS FOR N2:

T Vs H

           

T Vs S

             

T Vs Cp

                

PLOT FOR Co2:

T Vs H

               

T Vs S

                

T Vs Cp

             

CALCULATED MOLECULAR WEIGHT OF ALL SPECIES:

                                 

                                  

ERRORS FACED WHILE PROGRAMMING:

1.Forget to use the range operator in while calculating molecular weight.

     

2.Forget to use element wise operator while multiplying the formulas to calculate H,S and Cp

    

STEPS TAKEN TO OVERCOME THESE ERRORS:

1. By contacting the support engineer.

2. By looking at the red coloured popped hints at the command window.

CONCLUSION:

        Thus parsing of NASA thermodynamic data file to calculate the thermodynamic properties of various gas species and to save the plots is implemented successfully by writing the code logically in matlab.

REFERENCES:

https://in.mathworks.com/help/matlab/ref/cd.html

https://in.mathworks.com/matlabcentral/fileexchange/26139-molecular-weight-calculator-function

https://in.mathworks.com/help/matlab/ref/strtok.html?s_tid=srchtitle

 

 

 

 


Projects by Avinash A S

AIM:    To perform explicit dynamics simulation of a bullet penetrating into a bucket considering three different non-linear materials for the bucket and analyse their behaviour for the same bullet velocity. CASE SETUP:1 Assign the bucket material as Copp Read more

AIM:

    To perform an explicit dynamics simulation to simulate the machining operation with using a planer and observe their behaviour at 2 different velocities.

CASE SETUP:1Read more

AIM:    To perform the tension and torsion test on the specimen considering the case setups and find out the total deformation, equivalent stress and temperature of the specimen. TENSION TEST: A tensile test, also known as a tension test, is one Read more

AIM:    To perform transient structural analysis on a worm gear assembly. TRANSIENT STRUCTURAL ANALYSIS: Transient structural analysis is the process of calculating and determining the effects of loads and internal forces that are a function of time on a Read more

AIM:     To perform a transient structural analysis on a double universal joint with a spring using three different materials and then compare the results. TRANSIENT STRUCTURAL ANALYSIS: Transient structural analysis is the process of calculating and Read more

AIM:       To perform a transient structural analysis on a piston and cam mechanism model considering different case setups and then compare the results of the case setups. TRANSIENT STRUCTURAL ANALYSIS: Transient structural analysis is the process o Read more

AIM:      To perform structural analysisto simulate a sphere pressing on plate considering the case setup and experience the plastic deformation on the plate. STRUCTURAL ANALYSIS: Structural analysis is the process of calculating and determining the Read more

AIM:    To perform an structural analysis to simulate the bending of a wire considering three different case setups and compare their results. STRUCTURAL ANALYSIS: Structural analysis is the process of calculating and determining the effects of loads and Read more

AIM:     To perform structural analysis of weld joints for three different case setups and compare the directional deformation, equivalent stress and equivalent strain for all the cases STRUCTURAL ANALYSIS: Structural analysis is the process of calcu Read more

AIM:      To perform a static structural analysis on the Railwheel and Track setup considering the following case setups. STRUCTURAL ANALYSIS: Structural analysis is the process of calculating and determining the effects of loads and internal forces Read more


Loading...

The End