- To write a program in Matlab to parse the Nasa thermodynamic data file and Nasa polynomials.
- To extract coefficients of each species from the data file and calculate thermodynamic properties such as Specific heat , Enthaly & Entropy .
- To calculate the molecular weight of each species in the data .
- Automating the program in a way that user inputs only the required specie name and all the data is generated and stored for further evaluation.
- Nasa came up with polynomials that can be used to evaluate thermodynamic properties such as Cp,h and s.
- This data is also used for combustion calculations under combustion chemistry research.
- The Nasa Polynomials are :
- CpR=a1+a2â‹…T+a3â‹…T2+a4â‹…T3+a5â‹…T4
- HRT=a1+a2â‹…T2+a3â‹…T23+a4â‹…T34+a5â‹…T45+a6T
- SR=a1â‹…lnâ‹…T+a2â‹…T+a3â‹…T22+a4â‹…T33+a5â‹…T44+a7
- T = Temperature
- Cp = Specific Heat
- H = Enthalpy
- S = Entropy
- R = 8.314 Universal gas constant k−1⋅Mol−1
- a1, a2, a3, a4, a5, a6, and a7 are the numerical coefficients supplied in Nasa thermodynamic data file
- The Nasa thermodynamic data to parse is given here : http://combustion.berkeley.edu/gri_mech/version30/files30/thermo30.dat
- A snapshot of The first species of the data :
- Here in the first line the starting letter / letters given is the name of the species . In this case it is oxygen denoted by O . And towards the end of the first line after G , Local temperature ranges are given 200,3500 and 1000.
- Where 200k to 1000k is Local low temperature range and 1000k to 3500k is Local high temperature range.
- These Local temperature ranges vary for every species in a similar fashion explained above.
- From the second line a few numbers are given till the 4th line , they are 14 numerical coefficients we need to pass in the above polynomial equations in place of a1 ,a2 ,a3 .. likewise .
- The first 7 numbers starting on the second line of each species entry (five of the second line and two of the third line) are the seven coefficients (a1 through a7 , respectively) for the local high temperature range 1000k to 3500k .
- The following seven numbers are coefficients (a1 through a7 , respectively) for the Local low temperature range 200k to 1000k.
Matlab Project code : Here (Permission required) #
- Procedure of the Main program in matlab :
- With a fair understanding of the above terminologies of the data and polynomials , the procedure of the Matlab program is explained here .
- The Thermodynamic data file is opened using fopen command and read condition is assigned .
- Then the program uses input command to ask the user To enter the species name for which the data is to be generated .
-
- Then a while loop is created to read the entire file and from fgetl all the lines are stored in a cell variable tline.
- The lines are split to extract the first letter/letters of every specie line throughout the data file and the input command and the line where all the species name are stored is compared using strcmpi command (to accept capital and small case letters)
- Then after accessing the Specie line , the Local temperature range is extracted using a logic(explained ahead) and store it in variables .
- (data from matlab workspace for specie O2)
- After accessing the temperature ranges , the 3 lines ahead containing 14 coefficients is accessed using fgetl and converted into strings to apply the logic which is applied similarly for extracting the local temperature ranges .
- The logic is :
- Each Co-effcient numbers has the letter E . Now , using strfind command to get the location of the letter E , in the particular line , then a variable stores the coefficient number from start till the position of E and +3 . as all coefficients have 3 values after the letter E.
- similarly the next coefficient number extraction begins from The position E of the last number + 4th value , which is the beginning of the new coefficient number till the new position of E + 3 .
- The program part:
% Part of the code for explanation of logic
p = fgetl(fid); % getting the coefficient line
p = string(p); % converting to string
p1 = strfind(p(1),'E'); %finding the locations of E throughout the line
%output for p1 = [ 12 27 42 57 72]
% for 5 coefficients 5 E positions .
% getting the coefficient and converting it into numerical value of class
% double .
h1 = p{1}(1:p1(1)+3);
h1 = str2double(h1)
% Output of h1 = 2.9266 % for species N2 .​
- Similarly all 14 coefficients are extracted.
-
- (coefficients extracted in workspace of matlab for Co2)
- It proceeds within the counter created before while loop .
- The h variable refers to high temperature Coefficient and l variable refers to low temperature coefficient .
- Calculating Cp , h & s for Local low temperature range and Local high temperature range .
- The Temperature range extracted above are stored in variables called Local_low_temp, Local_mid_temp & Local_high_temp .
- They are now stored in another variable using linspace command for the Lower temperature range for eg : 200k - 1000k and for higher temp range 1000k - 3500k .
- The linspace command helps create a range of increasing values from 200k - 1000k and similarly for 1000k - 3500k . (The temp. ranges used here are just for example , actual ranges vary as per species)
- Now the Polynomials are programmed in matab and for one set of Polynomial equations Local low temperature range is passed . for eg : 200k - 1000k . Now these low temp range polynomials use last seven Coefficient numbers of every species as specified earlier
- And another set of Polynomial equations which pass on the higher Temp coefficients which are the first seven coefficients of every species. and local high temperature range 1000k-3500k is passed in this .
- Program :
-
Temp_lower = linspace(Local_low_temp,Local_mid_temp,500);
Temp_higher = linspace(Local_mid_temp,Local_high_temp,500);
T_l = Temp_lower;
cp_l = (l1 + l2.*T_l + l3.*T_l.^2 + l4.*T_l.^3 + l5.*T_l.^4).*R;
h_l = (l1 + l2.*(T_l)./2 + l3.*(T_l.^2)./3 + l4.*(T_l.^3)./4 + l5.*(T_l.^4)./5 + l6./T_l).*R.*T_l;
s_l = (l1.*log(T_l) + l2.*T_l + l3.*(T_l.^2)./2 + l4*(T_l.^3)./3 + l5*(T_l.^4)./4 + l7).*R;
T_h = Temp_higher;
cp_h = (h1 + h2.*T_h + h3.*T_h.^2 + h4.*T_h.^3 + h5.*T_h.^4).*R;
h_h = (h1 + h2.*(T_h)./2 + h3.*(T_h.^2)./3 + h4.*(T_h.^3)./4 + h5.*(T_h.^4)./5 + h6./T_h).*R.*T_h;
s_h = (h1.*log(T_h) + h2.*T_h + h3.*(T_h.^2)./2 + h4*(T_h.^3)./3 + h5*(T_h.^4)./4 + h7).*R;
- Each variable for cp , h and s stores 500 values .
- The Specific heat , Enthalpy and Entropy values are calculated for the species the user inputs .
- Plot Thermodynamic properties :
- The coefficients extracted from the data are passed into the polynomials , which then calculates the values for cp , h & s for local low temp range and local high temp range.
- The plot is created in two parts , one for low temperature range with low temperature coefficient and another for high temperature range with high temperature coefficients .
- The plots are then combined from the local temperature range (low temperature : high temperature) specific for each species and seperate colors are given to identify the low temp range with low temp coefficient and another high temp range with high temp coefficient.
- The plot is created for Specific heat vs Local Temperature for both ranges and coefficients . and similarly for Enthalpy vs Local temperature and Entropy vs Local temperature range.
- program :
-
plot(T_l,cp_l, 'linewidth',3,'color','r');
hold on
plot(T_h,cp_h, 'linewidth',3,'color','g');
xlabel('Local Temperature');
ylabel('Specific Heat');
title(sprintf(' Specific Heat vs Temperature for %s ', ask ));
legend('Low temp cp','High temp cp','location','northwest');
grid on
- Similarly for enthalpy and entropy different figures are plotted
- The output plots generated for O2 are :
- The plots are saved using saveas command in the directory set and are stored in the folder created for the specific species. this is generated automatically.
- Calculating Molecular weight of each species :
- The species contains the elements ['H' 'C' 'N' 'O' 'AR']
- Hence only those are used , and their atomic mass value is stored as an array [1 12 14 16 40]
- To calculate molecular weight of each species the logic used in shown in the code.
- It will calculate the molecular wt of the species the user inputs and store the value in a text file in the folder where the plots for Cp , h and S are saved . (ahead in the report)
- Program :
-
Species = ['H' 'C' 'N' 'O' 'AR'];
atomic_mass = [1 12 14 16 40];
mw = 0;
for i = 1:length(ask)
for j = 1:length(Species)
if strcmpi(ask(i),Species(j))
mw = mw + atomic_mass(j);
position = j;
end
end
n = str2double(ask(i));
if n > 1
mw = mw + atomic_mass(position)*(n-1);
end
end
fprintf('Molecular weight of %s is ', ask)
fprintf('%dn', mw)
fid1 = fopen('Molecular weight.txt','w');
fprintf(fid1,'Molecular weight of %s is %d', ask,mw);
fclose(fid1);
- The molecular weight along with saving in txt file , also displays in the command window of the specie name user inputs .
- Mkdir matlab command creates new folder .
- Since the program is coded in a way that user inputs the required specie to obtain the data. The command mkdir is assigned the name of the species .
- Hence it creates the folders with the specie names whenever the user inputs and runs the program in the directory the main program and thermodynamic data file is stored .
- After creating the plots and getting molecular weight , the saveas command to save the plots and writing the molecular weight in the txt file and saving it in the same folder is done together at the end after specifying the directory and the name of the folders mkdir will create .
- At the end the directory is changed again to the original directory to run the program again with just one user input and automatically evaluate and theres no need to manually go back .
Folders containing the thermodynamic data for all species
All the folders are created in this directory and plots and molecular wt txt file is saved in there. Any specie the user inputs the folder is created and data is generated and stored.
- Folder containing plots & weight file for o2
- Folder containing plots & weight file for N2
- Plots for N2 and Co2 Species :
Plot for N2
Plot for Co2
- For more references regarding the Nasa thermodynamic data and polynomials :
- http://combustion.berkeley.edu/gri_mech/version30/files30/thermo30.dat
- http://combustion.berkeley.edu/gri_mech/data/nasa_plnm.html
- https://www.researchgate.net/publication/289281167_Thermochemical_Data_for_Combustion_Calculations
- https://www.grc.nasa.gov/WWW/CEAWeb/TP-2002-211556.pdf
- https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19720021272.pdf
THE END
keywords - MATLAB, CFD, MATLAB-BASICS, NUMERICAL-ANALYSIS