## 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.

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)

A VIEW OF NASA THERMODYNAMIC FILE:

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.

PLOTS FOR o2:

T Vs H

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

### CAR CRASH SIMULATION USING ANSYS WORKBENCH Avinash A S · 2020-02-27 12:32:21

AIM:      To perform a car crash analysis using ansys workbench. PROBLEM STATEMENT: To perform a parametric study using 3 different values of thickness for the car body. Here the car body remains stationary but the wall displaces to a distance of 500 Read more

### COPPER ROLLING SIMULATION IN ANSYS WORKBENCH Avinash A S · 2020-02-25 12:15:10

AIM:     To perform structural analysis to simulate the rolling operation on a copper workpiece of length 90mm. STRUCTURAL ANALYSIS: Structural analysis is the process of calculating and determining the effects of loads and internal forces on a struc Read more

### EXPLICIT DYNAMICS SIMULATION OF A BULLET PENETRATING INTO A BUCKET Avinash A S · 2020-02-12 11:33:34

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

### EXPLICIT ANALYSIS OF THE MACHINING WITH PLANER Avinash A S · 2020-02-09 05:39:33

AIM:

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

### TENSION AND TORSION TEST ON STEEL SPECIMEN USING ANSYS WORK BENCH Avinash A S · 2020-02-08 07:01:58

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

### TRANSIENT STRUCTURAL ANALYSIS ON WORM GEAR ASSEMBLY USING ANSYS WORKBENCH Avinash A S · 2020-02-07 13:06:11

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

### TRANSIENT STRUCTURAL ANALYSIS OF DOUBLE UNIVERSAL JOINT USING ANSYS WORKBENCH Avinash A S · 2020-02-07 11:24:22

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

### TRANSIENT STRUCTURAL ANALYSIS OF LONG PISTON WITH CAM USING ANSYS WORKBENCH Avinash A S · 2020-02-06 10:49:41

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

### ANALYSIS OF SPHERE PRESSING ON PLATE USING ANSYS WORKBENCH Avinash A S · 2020-01-29 08:41:21

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

### BENDING OF WIRE USING ANSYS WORKBENCH Avinash A S · 2020-01-28 09:18:24

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