Project 1 - Parsing NASA thermodynamic data

AIM: Write a MATLAB program that extracts the 14 coefficients and calculates the enthalpy, entropy, and specific heat for all the species in the data file (NASA Thermodynamic data).

 

Introduction:

NASA introduced polynomials that can be used to calculate the thermodynamic properties such as Specific heat(Cp), Enthalpy (H), and Entropy(S) of the elements and their species. They also have documented some coefficients that are required to evaluate these polynomials for 1000+ species of elements.

In order, to calculate the thermodynamic properties of such elements and their species, we will require the data file of NASA. Similarly, we will be using the technique of "file parsing" to read and evaluate the data.

 

What is file parsing?

Parsing in computer languages refers to the syntactic analysis of the input data normally from a file into its components parts to facilitate the writing of compilers and interpreters. In simple words, it refers to how information can be read from a text file or any type of file. And how that information can be used in a computer program.
By using file parsing, we will be able to read data from external files and we can include data in our calculation.
When we talk about data processing there is an interface between the program and a large amount of data. So how we can carry out file parsing?
We can use various efficient software or we can program it in a programming interface platform. As far as our project is concerned, we will be parsing technique in MATLAB i.e.(Programming interface platform).
 

Theory:

In our project, we will be calculating various thermodynamic properties of the elements and their species. So let's get a fundamental recap of such thermodynamic properties.

1. Specific Heat(Cp): Specific heat is defined as the amount of heat that must be added to one unit mass of the substance to cause an increase of one unit in temperature. The unit of specific heat is J/Kg K.

2. Enthalpy(H): Enthalpy is the sum of internal energy(U) and product of pressure and volume(PV). In simple terms, it tells how much heat and work was added or removed from the substance. It is similar to energy, but not the same. The unit of enthalpy is KJ/kg.

3. Entropy(S): Entropy is the measure of a system's thermodynamic 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 molecular disorder or randomness, of a system. The unit of entropy is KJ/Kg K.

 

Objectives: 

1. Write a MATLAB program to calculate the specific heat, enthalpy, and entropy of the species.

2. Calculate the molecular weight of each species and display it in the command window.

3. Plot the specific heat, enthalpy, and entropy for the local temperature range.

4.  Create separate folders for each species and save the plots as images.

 

Commands used in this project:

i] fopen - This MATLAB command opens the file or obtains information about open files.

ii] fgetl - This command reads the line from the file or removes newline characters.

iii] strfind - It searches for strings within other strings.

iv] str2double - It converts strings to double-precision values.

v] strsplit - This MATLAB command splits the string or character vector at a specified delimiter.

vi] linspace - It generates a linearly spaced vector.

vii] for - It repeats the loop for a specified number of times.

viii] if - It executes the statement "if" condition is true.

ix] sprintf - This MATLAB command formats data into string or character vector.

x] mkdir - This command makes a new folder.

xi] cd - It changes the current folder.

xii] saves - This MATLAB command saves figures to the specific file format.

Now we creating a function file for,' cp', 'H','s' 

function file with '.m' is created where we are going to write the formula for,' specific heat ', 'entropy' and enthalpy, where (a1,a2,.....a7) are coefficients of local high temperature where other values of (a8,a9,.....a14) are coefficients of local low temperature.

now we are going to create a 'for loop' for local temperature species. with 'for loop' we are going to create an 'if loop' inside a 'for loop'. the need of 'if loop' states the condition that if Temperature(T) exceeds the local temperature then "Cp, H, S'' will be calculated for local high temperature with their respective coefficients. while the else condition states that if the temperature doesn't exceed the local temperature(T) then the calculated values for the species will be for local low temperature with their dedicated local coefficients.

Formula :

MATLAB CODE:

% MATLAB CODE: for "Cp, H, S "

function [Cp,H,S] = file_func(T,R,local_mid_temp,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14)

for i = 1:length(T)
if T(i)>local_mid_temp
    Cp(i) = (a1 + a2*T(i) + a3*T(i)^2 + a4*T(i)^3 + a5*T(i)^4).*R;
    H(i) = (a1 + (a2*T(i))/2 + (a3*T(i)^2)/3 + (a4*T(i)^3)/4 + (a5*T(i)^4)/5 + (a6/T(i))).*T(i)*R;
    S(i) = (a1*log(T(i)) + a2*T(i) + (a3*T(i)^2)/2 + (a4*T(i)^3)/3 + (a5*T(i)^4)/4 + a7).*R;
   
else
    Cp(i) = (a8 + a9*T(i) + a10*T(i)^2 + a11*T(i)^3 + a12*T(i)^4).*R;
    H(i) = (a8 + a9*T(i)/2 + (a10*T(i)^2)/3 + (a11*T(i)^3)/4 + (a12*T(i)^4)/5 + (a13/T(i))).*T(i)*R;
    S(i) = (a8*log(T(i)) + (a9*T(i)^1) + (a10*T(i)^2)/2 + (a11*T(i)^3)/3 + (a12*T(i)^4)/4 + a14).*R;
end

end
end 

Function coding to get the molecular weight: 

i) A separate ".m file" for molecular weight calculations and a function is created.

ii) The inputs are given for atomic weight and species name.

iii) The values of elements are given. For later calculations of the molecular weight of species, these values are useful.

iv) After that, a "for loop" is created for calculating the molecular weight of the species individually.

Note: Two loops are created because the species are part of elements. The atomic weight of elements is a datum for the species, by which they are getting evaluated.  

The "strcmp" command is used to compare species with the elements and to calculate their molecular weight.

The "str2double" command is used to convert the strings into numbers

 MATLAB CODE MOLECULAR WEIGHT:

%MATLAB code: for molecular weight


function [molecular_weight_0] = molecular_func(species_name)
species = species_name;
elements = ['H','C','O','N','A','S'];
atomic_weight = [1.00784,12.0107,15.999,14.0067,39.948,32.065];
molecular_weight_0 = 0;
for i = 1:length(species)
    for j = 1:length(elements)
        if strcmp(species(i),elements(j))
           
   molecular_weight_0 = molecular_weight_0 + atomic_weight(j);
   x = j;
   
        end
    end
    n = str2double(species(i));
    if n>1
       molecular_weight_0 = molecular_weight_0 + (atomic_weight(x)*(n-1));
    end
end

 Main Program:

In this program, we will code all the step-by-step syntax to execute the problem statement with appropriate commands and logic.

Step by step procedure:

a) First we will open the "THERMO.dat" data file by using the "fopen" command 

b) Then the "fgetl" command used gets the first line of the file(z). Similarly, we will get the second line for the global temperature header line, and assign a variable "temp_header" to it.

c) From this temp_header we are going to get the string of global low, mid and high temperature by using "strfind" command and we will search for '.' (Fullstop) position in that respective line. "K" is the variable assigned to find the string.

d) After getting the position of "." in temp_header, we will be calculating global low temperature by inserting the position of the "." (full stop) and manipulating the position to get the values before and after the "." (full-stop). Similarly, we will calculate global mid and high temperature and will assign an independent variable for the calculation. Finally, converting the string into a number by using the "str2double" command. Now, the global temperature value is found out.

e) We will create a "for loop" by inserting the "fgetl" command and iterate it from the 3rd to the 5th line To skip unwanted lines.

f) As stated in the data file, we have 53 species, therefore assigning the "num_species = 53".

"For loop started"

g) create a loop that will be having the syntax that will execute the calculations for every species and will be closed at the end of the program.

h) The 6th line, is the header line of the species. Again using the "fgetl" command to get the header line printed and assigning a variable to it. Based on the logic for the global temperature calculations i.e.[step (d)], calculating the local low, mid and high temperature.

i) To get, the species name printed for all the iterations, splitting the string using the "strsplit" command and assign an independent variable to it. Then string gets split into cell arrays from which have to extract the first cell which shows the name of the species.

j) calculated the coefficients of species that are equally spaced in 3 consecutive lines. For this, getting the first line of coefficients by using the "fgetl" command and assign a variable to the syntax. Accordingly, find the common character in every coefficient, which nothing but "E".

k) found the position of "E" by using the "strfind" command in the first coefficient line and after getting that modify the values such that the whole coefficient should be found. When the extracted coefficient is in the form of the string, then convert the string into a number. using the "strdouble" command for the conversion. Similarly, executed the syntax for each coefficient and the other two coefficient lines.

l) Providing all the coefficient calculations in our program. Then calculate the specific heat, enthalpy, and entropy for the given i.e. local temperature range. So on define the local temperature range using the "linspace" command and inserting universal gas constant. Then on calling the function of Cp, H, and S, which have been created above. Verifying our function working by pressing the "F5" key. Hundreds of different values of Cp, H, and S  for hundred linearly spaced different local temperature range.

m) Before plotting the values of Cp, H, and S, create a directory (Folder) in which all the images of the plots will be saved accordingly. Here, on the use of the "mkdir" and "cd.." command.

n) plotting the values of Cp, H, and S w.r.t temperature variation, label the x and y-axis, and save the images of the plot by using the "saveas" command.

o) calculated the molecular weight by using the molecular weight function, which has created previously as an independent ".m file". then the same function will be called and will print the molecular weight by using the "sprintf" command with %s %d by which the value of molecular mass will be displayed.

NOTE: The for loop runs in a way that if there is a numeric value in the species name, then it will extract that value, converting it in number and then adds the previous character preceding the number, that number of times to get the molecular weight.

E.g: For O2, it will extract 2 from the string, convert it into the numeric value and add O two times (O+O), which will result in the molecular weight of O2.

"For loop ended"

p) Press the "F5" key to run the program and after execution, open the directory stated in the program to visualize the results. 

 MATLAB CODE(MAIN)

clear all
close all
clc

z = fopen('THERMO.dat','r');
fgetl(z);
temp_header = fgetl(z);
K = strfind(temp_header,'.');

% To calculate global low temperature.
% e = temp_header((7-3):7+3)
e = temp_header((7-3):K(1)+3);
global_low_temp = str2double(e);

% To calculate global high temperature.
% f = temp_header(10:17+3)
f = temp_header(K(1)+6:K(2)+3);
global_mid_temp = str2double(f);

% To calculate global mid temperature.
% g = temp_header(20+3:(27+3))
g = temp_header(K(2)+6:K(3)+3);
global_high_temp = str2double(g);

% For loop for neglecting unwanted lines.
for i = 3:5
    skip_lines = fgetl(z) ;
end

% Total number of species.
num_species = 53;

for j =1:num_species;
% To get the species and local low, mid and high temperature.
species_header = fgetl(z);
A = strfind(species_header,'.');
% To get the local low temperature
% local_temp_1 = species_header(52-3:52+3)
local_temp_1 = species_header(A(1)-3:A(1)+3);
local_low_temp = str2double(local_temp_1);

% To get the local high temperature
% local_temp_2 = species_header(52+6:62+3)
local_temp_2 = species_header(A(1)+6:A(2)+3);
local_high_temp = str2double(local_temp_2);

% To get the local mid temperature
% local_temp_3 = species_header(62+6:72+3)
local_temp_3 = species_header(A(2)+6:A(3)+3);
local_mid_temp = str2double(local_temp_3);

% To get the species name
q = strsplit(species_header,' ');
species_name = q{1};

% To calculate all coefficients in line 1.
coeff_line_1 = fgetl(z);
B = strfind(coeff_line_1, 'E');

% Coefficients from line 1
% 1st coefficient
a1 = coeff_line_1(1:B(1)+3);
a1 = str2double(a1);
% 2nd coefficient
a2 = coeff_line_1(B(1)+4:B(2)+3);
a2 = str2double(a2);
% 3rd coefficient
a3 = coeff_line_1(B(2)+4:B(3)+3);
a3 = str2double(a3);
% 4th coefficient
a4 = coeff_line_1(B(3)+4:B(4)+3);
a4 = str2double(a4);
% 5th coefficient
a5 = coeff_line_1(B(4)+4:B(5)+3);
a5 = str2double(a5);

% To calculate all coefficients from line 2.
coeff_line_2 = fgetl(z);
C = strfind(coeff_line_2, 'E');

% Coefficients from line 2
% 6th coefficient
a6 = coeff_line_2(1:C(1)+3);
a6 = str2double(a6);
% 7th coefficient
a7 = coeff_line_2(C(1)+4:C(2)+3);
a7 = str2double(a7);
% 8th coefficient
a8 = coeff_line_2(C(2)+4:C(3)+3);
a8 = str2double(a8);
% 9th coefficient
a9 = coeff_line_2(C(3)+4:C(4)+3);
a9 = str2double(a9);
% 10th coefficient
a10 = coeff_line_2(C(4)+4:C(5)+3);
a10 = str2double(a10);

% To calculate all coefficients from line 3.
coeff_line_3 = fgetl(z);
D = strfind(coeff_line_3, 'E');

% Coefficients from line 3
% 11th coefficient
a11 = coeff_line_3(1:D(1)+3);
a11 = str2double(a11);
% 12th coefficient
a12 = coeff_line_3(D(1)+4:D(2)+3);
a12 = str2double(a12);
% 13th coefficient
a13 = coeff_line_3(D(2)+4:D(3)+3);
a13 = str2double(a13);
% 14th coefficient
a14 = coeff_line_3(D(3)+4:D(4)+3);
a14 = str2double(a14);


% To calculate specific heat, Enthalpy and Entropy
% Local temperature range
T = linspace(local_low_temp, local_high_temp,100);
% Universal gas constant
R = 8.314;
% Function to calculate specific heat, enthalpy and entropy.
[Cp,H,S] = file_func(T,R,local_mid_temp,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14);

% Creating folder
mkdir([species_name])
cd([species_name])

% Plotting
figure(1)
plot(T,Cp,'linewidth',3);
xlabel('Temperature [K]');
ylabel('Specific Heat [KJ/Kg K]');
grid on
title('Temperature vs Specific heat')
saveas(1,'Specific Heat','jpg');

figure(2)
plot(T,H,'linewidth',3);
xlabel('Temperature [K]');
ylabel('Enthalpy [KJ/Kg]');
grid on
title('Temperature vs Enthalpy')
saveas(2,'Enthalpy','jpg');

figure(3)
plot(T,S,'linewidth',3);
xlabel('Temperature [K]');
ylabel('Entropy [KJ/Kg K]');
grid on
title('Temperature vs Entropy')
saveas(3,'Entropy','jpg');
cd ..
% To calculate the molecular weight
[molecular_weight] = molecular_func(species_name);
sprintf('Molecular weight of %s is %d', species_name,molecular_weight)
end

fclose(z);

 All the Cp, h, and S, plots are stored in their species folder

 Results:

1) plots for O2, N2 andCO2 species

A. plots for O2

 B. plots for N2

  C. plots for CO2

WorkSpace ss

Error and correction

here, the error occurred is by function name as the function file is named as file_func1 but while creating the function it is named as file_func which is not recognized by Matlab so by changing the file name from file_func1 to file_func then the error is corrected. and after running this it shows zero error.

Conclusion

In this project, calculation of thermodynamic properties of the respective species is calculated by Cp, H, and S. The logic for approaching the solution is that the problem is broke down into simple intervals and code accordingly, it is good for the ease to approach the whole code for single species and make slight changes by adding for loops and execute the program. cracking the data file with sequence or line by line, it is seen that having MATLAB basic knowledge the execution of the project can be done greatly. Therefore, the given thermodynamic properties with respective plots and molecular weight calculations were resulted by using MATLAB programming.

References

1) https://uk.mathworks.com/matlabcentral/answers/7262-best-way-to-parse-data-from-a-large-mixed-format-text-file

2) https://skill-lync.com/projects/FILE-PARSING-WITH-MATLAB-93659

3) https://www.govinfo.gov/content/pkg/GOVPUB-C13-8977bb85b16c045950bc566ad4039657/pdf/GOVPUB-C13-8977bb85b16c045950bc566ad4039657.pdf

4) https://nvlpubs.nist.gov/nistpubs/Legacy/TN/nbstechnicalnote129.pdf

5)https://uk.mathworks.com/matlabcentral/answers/123314-how-to-read-numeric-data-from-text-file-and-plot-graph-for-given-numeric-data

6)https://pubs.acs.org/doi/10.1021/acsomega.9b02488#:~:text=These%20results%20indicate%20that%20our,K%E2%80%931%20at%20298%20K.

7)https://www.researchgate.net/publication/283693879_Thermodynamic_and_Transport_Properties_of_N2O2_Mixtures_at_Different_Admixture_Ratios