Project 1 - Parsing NASA thermodynamic data
- Aim:
To write a code in Matlab to parse the NASA thermodynamic data file and then calculate the thermodynamic properties of various gas species.
- Objective:
1. To write a function that extracts the 14 coefficients and calculates the enthalpy, entropy, and specific heats for all the species in the NASA thermodynamic data file "Thermo.dat
The NASA polynomials have the form as shown below:
Cp/R = a1+ a2 T + a3 T^2 + a4 T^3+ a5 T^4
H/RT= a1+ a2 T /2 + a3 T 2 /3+ a4 T3 /4 + a5 T4 /5+ a6/T
S/R =a1 lnT+ a2 T + a3 T^2 /2 + a4 T3 /3+ a5 T4 /4 + a7
where R is the universal gas constant, T is temperature
2. To calculate the molecular weight of each species and display it in the command window.
3 To plot the Cp. Enthalpy, and Entropy for the local temperature range (low temperature: high temperature) specific for each species.
4. To save the plots as images with appropriate names and dump them in separate folders for each species with suitable name automatically
- Introduction:
File parsing:
File parsing refers to how information can be read from a particular text file or any type of file and how that Information can be used in a computer program
Step 1-Define the functionStep 2-Create on inputParser objectStep 3- Add inputs to the scheme Step 4- Set properties to adjust parsing Step 5-Parse the inputs. Step 6-Use the inputs in function Step 7-Calling function
NASA Thermodynamic File
THERMO dat file is a data file containing 14 coefficients for 53 different elements given by NASA Those coefficients are utilized to calculate the specific heat, enthalpy & entropy for a particular element. Each element has its own range of local temperature. Below screenshot of the THERMO dat file is displayed.
List of 53 Elements contained in THERMO.dat file:
C C2H AR C2H2 C2H3 C2H4 CH C3H7 C3H8 C2H5 C2H6 CH2 CH20 CH2OH CH3 CH2(s) CH2CHO CH2CO CH4 CN CO CH3CHO CH30 CH3OH H2CN H20 H202 H2 CO2 HHCNO HCO HCNN HCCO HCN CCOR HOCN N N2 HNCO HNO HO2 NH3 NH2 NNH NH N20 NCO 02 OH NO NO2
Enthalpy, Entropy, and Specific Heat:
Enthalpy, entropy, and specific heat, of various species values at different temperatures are required extensively to analyze the thermodynamic process
For example, if the combustion process is considered during this, these quantities are required to observe how the combustion is going on and how efficient it will be. It eventually allows seeking the scope of improvement in the process. Similarly in the other processes like refrigeration, air conditioning, chemical kinetics, nozzle flow, detonation process, etc, these quantities are widely involved to observe the phenomenon.
NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H, and 5 They have documented the coefficients that are required to evaluate these polynomials for 1000+ species
The NASA polynomials have the form
Cp/R = a1+ a2 T + a3 T^2 + a4 T^3+ a5 T^4
H/RT =a1+ a2 T /2 + a3 T 2 /3+ a4 T3 /4 + a5 T4 /5+ a6/T
S/R =a1 lnT+ a2 T + a3 T^2 /2 + a4 T3 /3+ a5 T4 /4 + a7
Where R is the universal gas constant, T is the temperature and al, a2 a3, 04, 05, a6, and a7 are the numerical coefficients supplied in the thermodynamic file. The FIRST 7 coefficients are HIGH-temperature coefficients (above 1000 K) and the SECOND 7 coefficients are LOW-temperature coefficients (below 1000 K).
- Coding Functions:
1.Specific Heat:
It is the amount of energy that must be added, in the form of heat, to one unit of mass of the substance in order to cause an increase of one unit in its temperature
Cp = (a1+ a2 T + a3 T^2 + a4 T^3+ a5 T^4)*R
function Cp=specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,R,global_medium_temp)
%compare with global medium temp
%if temp is less than or equal to global medium temp, minimum value of coefficient are taken
if temp<=global_medium_temp
Cp=R.*(a8+(a9*temp)+(a10*temp.^2)+(a11*temp.^3)+(a12*temp.^4));
else
Cp=R.*(a1+(a2*temp)+(a3*temp.^2)+(a4*temp.^3)+(a5*temp.^4));
end
2. Entropy:
Entropy, the measure of a system's thermal energy per unit temperature that is unavailable for doing useful work Because work is obtained from ordered molecular motion, the amount of entropy is also a measure of the molecular disorder, or randomness, of a system. The concept of entropy provides deep insight into the direction of spontaneous change for many everyday phenomena.
S =(a1*lnT+ a2 T + a3 T^2 /2 + a4 T3 /3+ a5 T4 /4 + a7)*R
function S=entropy_calc(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,R,global_medium_temp)
%compare with global medium temp
%if temp is less or equal to global medium temp minimum value
if temp <= global_medium_temp
S= R.*((a8.*log(temp))+(a9.*temp)+(a10*(temp.^2)/2)+(a11*(temp.^3)/3)+(a12*(temp.^4)/4)+a14);
else
S=R.*((a1.*log(temp))+(a2.*temp)+(a3*(temp.^2)/2)+(a4*(temp.^3)/3)+(a5*(temp.^4)/4)+a7);
end3. Enthalpy:
Enthalpy, a property of a thermodynamic system, is equal to the system's internal energy plus the product of its pressure and volume. In a system enclosed so as to prevent the mass transfer, for processes at constant pressure the heat absorbed or released equals the change in enthalpy.
H =(a1+ a2 T /2 + a3 T 2 /3+ a4 T3 /4 + a5 T4 /5+ a6/T)*R*T
function H=enthalpy_calc(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,temp,R,global_medium_temp)
%compare with global medium temp
%if temp is less or equal to global medium temp minimum values are taken
if temp <= global_medium_temp
H=R.*temp.*(a8+(a9.*(temp/2))+(a10*(temp.^2)/3)+(a11*(temp.^3)/4)+(a12*(temp.^4)/5)+(a13./temp));
c
else
H=R*temp.*(a1+(a2.*(temp/2))+(a3*(temp.^2)/3)+(a4*(temp.^3)/4)+(a5*(temp.^4)/5)+(a6./temp));
end
4. Molecular Weight:
The molecular mass is the mass of a given molecule: it is measured in daltons. Different molecules of the same compound may have different molecular masses because they contain different isotopes of an element
Function for Molecular Weight,
function Molecularweight=Molecular_weight(Species_SP)
%Standard Atomic Name
Atomic_name = ['H','C','O','N','A','S'];
Atomic_weight=[1.008 12.011 15.999 14.007 39.948 32.00];
Molecularweight=0;
for i=1:length(Species_SP)
for j=1:length(Atomic_name)
if strcmp(Species_SP(i), Atomic_name(j))
Molecularweight=Molecularweight + Atomic_weight(j);
position=j;
end
end
n=str2num(Species_SP(i));
if n>1
Molecularweight=Molecularweight + Atomic_weight(position)*(n-1);
end
end
fprintf('molecular weight %s : ', Species_SP);
fprintf('%f',Molecularweight);
end- Program:
clear all
close all
clc
% Reading Data file- NASA Thermodynamic
% Open a file & 'r' Open the file for reading
f1=fopen('THERMO.dat','r');
gt=fgetl(f1);
%Obtaing Temperatures-by Reading Header Information
global_temp=fgetl(f1);
% Spliting
t=strsplit(global_temp);
% Converting string cells to numbers using str2num,
global_low_temp=str2num(t{2});
global_medium_temp=str2num(t{3});
global_high_temp=str2num(t{4});
% Temperature range & Gas Constant
%Temperature Range assumed
temp=linspace(global_low_temp,global_high_temp, 1000);
%Universal Gas Constant J/(mol-K)
R = 8.314;
% skipping three comment lines in the file
for g = 1:3
gt=fgetl(f1);
end
for h = 1:53
gt=fgetl(f1);
A=strsplit(gt,' ');
Species_SP=(A{1});
posdot=strfind(gt,'.');
lowtemp = gt(posdot(1)-3:posdot(1)+3);
lowtemp =str2num(lowtemp);
line_1=fgetl(f1);
a=strfind(line_1, 'E');
% Extracting higher and lower temperature coefficients
%%Extracting FIRST 7 coefficients i.e HIGH-temperature
a1=line_1(1:a(1)+3);
a1=str2num(a1);
a2=line_1(a(1)+4:a(2)+3);
a2=str2num(a2);
a3=line_1(a(2)+4:a(3)+3);
a3=str2num(a3);
a4=line_1(a(3)+4:a(4)+3);
a4=str2num(a4);
a5=line_1(a(4)+4:a(5)+3);
a5=str2num(a5);
line_2=fgetl(f1);
b=strfind(line_1, 'E');
a6=line_2(1:b(1)+3);
a6=str2num(a6);
a7=line_2(a(1)+4:b(2)+3);
a7=str2num(a7);
% Extracting SECOND 7 coefficients LOW-temperatur
a8=line_2(b(2)+4:b(3)+3);
a8=str2num(a8);
a9=line_2(b(3)+4:b(4)+3);
a9=str2num(a9);
a10=line_2(b(4)+4:b(5)+3);
a10=str2num(a10);
line_3=fgetl(f1);
c=strfind(line_1, 'E');
a11=line_3(1:c(1)+3);
a11=str2num(a11);
a12=line_3(c(1)+4:c(2)+3);
a12=str2num(a12);
a13=line_3(c(2)+4:c(3)+3);
al3=str2num(a13);
a14=line_3(c(3)+4:c(4)+3);
a14=str2num(a14);
% Specific heat calculation
Cp=specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,R,global_medium_temp);
% Enthalpy calculation
H=enthalpy_calc(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,temp,R,global_medium_temp);
% Entropy calculation
S=entropy_calc(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,R,global_medium_temp);
% Calculating molecular weight of species
Molecularweight= Molecular_weight(Species_SP);
% Read & Write Molecular list
Molecular_list=fopen('Molecular weight.txt', 'a+');
fprintf( Molecular_list,'Molecular weight of %s is %d ',Species_SP,Molecularweight);
fclose(Molecular_list);
% plotting (CP,H,S,Temp, Species_SP)
%Saving as a folder
cur_dir=pwd;
mkdir(Species_SP);
cd(Species_SP);
%Plot 1, Specific heat (CP) vs Temperature Range (Temp)
figure(1)
plot(temp,Cp, 'linewidth',2,'color','r');
xlabel('Temperature [K]');
ylabel('Specific_Heat [KJ]');
title(sprintf('Varying Specific Heat with Temperature Range of %s',Species_SP))
grid on
saveas(figure(1), 'Specific Heat.jpg');
%Plot 2, Enthalpy (H) vs Temperature Range (Temp)
figure(2)
plot(temp,H,'linewidth',2,'color','g');
xlabel(' Temperature [K]');
ylabel('Enthalpy in [k3/(mol)]');
title(sprintf( 'Varying Enthalpy with Temperature Range of %s', Species_SP));
grid on;
saveas(figure(2), 'Ethalpy.jpg')
%Plot 3, Entropy (5) vs Temperature Range (Temp)
figure(3)
plot(temp,S,'linewidth',2,'color','b');
xlabel('Temperature[K]');
ylabel('Entropy in (K3/(mol-)');
title(sprintf('Varying Entropy with Temperature Range %s', Species_SP));
grid on;
saveas (figure(3), 'Entropy.jpg');
cd(cur_dir);
end- Explaination:
1. start with clear all, close all, clc
2.reading the text file "Thermo.dat' file and opening it by using 'fopen' and 'r' as reading command
3. For reading the line in a text file using 'fgetl'.
4. To obtain the temperature range from the Thermo dat file, using 'strsplit' to split the temperature into a particular string cell.
5. Used 'str2num' for converting string cell to numbers.
6. Divided the temperature into a string and then converted to a number so that temperature values can be used separately to solve the NASA thermodynamics.
7. Using for & end' loops to skip the next three comment lines in the text file.
8. Reading all the 53 species to plot specific heat(cp), enthalpy(H), and entropy(S), with temperature values, and to find the molecular weights
9. Used 'for' loop to read 53 species.
10. Used 'strfind' to find string
11. Using 'linspace' specified the temperature range from low to high
12. Adding universal gas constant value 'R'
13. Extracted higher and lower temperature coefficient by reading line by line from the species
14. Colling specific heat, enthalpy, entropy function codes
15. Calculating Molecular weight.
16. With the help pwd mkdir and cd command we create folder for every element.
17. Plot graph specific value vs Temprature.
18. Plot graph Enthalpy vs Temperature.
19. Plot graph Entropy vs Temperature.
Syntax Used during Coding
1. Fopen: Open a file or obtain information about open files
2. r: Open the file for reading
3. fgetl: Read line from file
4. strsplit: Split string or character vector at specified delimiter
5. str2num: Convert character array or string to numeric array
6. strfind: Find one string within another
7. linspace: Generate linearly spaced vectors
8. fprintf: Write formatted data to file
9. pwd: Identify current directory
10. mkdir: Make new directory
11. cd: Change working directory
- Output:
Molecular weight AR: 39.9480
Molecular weight C: 12.0110
Molecular weight C2H: 25.0300
Molecular weight C2H2: 26.0380
Molecular weight C2H3: 27.0460
Molecular weight C2H4: 28.0540
Molecular weight C2H5: 29.0620
Molecular weight C2H6: 30.0700
Molecular weight C3H7: 43.0890
Molecular weight C3H8: 44.8970
Molecular weight CH: 13.0198
Molecular weight CH2: 14.0270
Molecular weight CH2CHO: 43.0450
Molecular weight CH2CO: 42.0370
Molecular weight CH20: 30.0260
Molecular weight CH2OH: 31.0340
Molecular weight CH3: 15.0350
Molecular weight CH3CHO: 44.0530
Molecular weight CH30: 31.0340
Molecular weight CH3OH: 32.0420
Molecular weight CH4: 16.0430
Molecular weight CN: 26.0180
Molecular weight CO: 28.0100
Molecular weight CO2: 44.0090
Molecular weight H: 1.0080
Molecular weight H2: 2.0160
Molecular weight H2CN: 28.0340
Molecular weight H20: 18.0150
Molecular weight H202: 34.0140
Molecular weight HCCO: 41.0290
Molecular weight HCCOH: 42.0370
Molecular weight HCN: 27.0260
Molecular weight HCNN: 41.0330
Molecular weight HCNO: 43.0250
Molecular weight HCO: 29.0180
Molecular weight HNCO: 43.0250
Molecular weight HNO: 31.0140
Molecular weight H02: 33.0060
Molecular weight HOCN:43.0250
Molecular weight N: 14.0070
Molecular weight N2: 28.0140
Molecular weight N20: 44.9130
Molecular weight NCO: 42.0170
Molecular weight H: 15.0158
Molecular weight NH2: 16.0230
Molecular weight NH3: 17.0310
Molecular weight NNH: 29.0220
Molecular weight NO: 30.0060
Molecular weight NO2: 46.0050
Molecular weight O: 15.9998
Molecular weight O2: 31.9980
Molecular weight OH: 17.0070
Cp vs Temp, H vs Temp, S vs Temp graph of Some Elements:
For O:
For Ar:
For H2:
- Conclusion:
Using MATLAB code, the necessary coefficients required to evaluate various thermodynamic parameters for different species are extracted from NASA thermodynamic data file. The specific heat, enthalpy, and entropy for all the species are calculated and plotted against the local temperature range The molecular weight is calculated for all the species.
google drive link given below for folder,
https://drive.google.com/drive/folders/1OcjLKWHx9Gzly9bBknARAlrWfO9j346t?usp=sharing
- Reference:
https://in.mathworks.com/help/matlab/ref/mkdir.html
https://in.mathworks.com/matlabcentral/answers/317006-create-new-folders-by-mkdir
https://in.mathworks.com/help/matlab/ref/saveas.html
https://in.mathworks.com/help/matlab/ref/fprintf.html
https://in.mathworks.com/help/matlab/ref/sprintt.html
https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19940013151.pdf