BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

Project 1 - Parsing NASA thermodynamic data

RRahul Bhagat6 Feb 2023

Aim: 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.

Objective:

Writing a function that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the data file. By reading this file NASA Thermodynamic data file. (Tip: use fopen and fgetl). The formulae are shown below.

Here R is the universal gas constant, T is the temperature

The FIRST 7 coefficients are HIGH-temperature coefficients and the SECOND 7 coefficients are LOW-temperature coefficients.

By using LOCAL TEMPERATURE RANGES which is unique for each individual species, not the GLOBAL TEMPERATURE RANGE.

  1. Calculate the molecular weight of each species and display it in the command window.
  2. Plotting the Cp, Enthalpy and Entropy for the local temperature range (low temp high temp) specific for each species.
  3. Saving the plots as images with appropriate names and dump them in separate folders for each species with suitable name.

Theroy:

Specific heat:

The specific heat capacity of a substance is the heat capacity of a sample of the substance divided by the mass of the sample. Informally, 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. The SI unit of specific heat is joule per kelvin and kilogram, J/.

The specific heat constants for constant pressure and constant volume processes are related to the gas constant for a given gas. This rather remarkable result has been derived from thermodynamic relations, which are based on observations of physical systems and processes. Using the kinetic theory of gases, this same result can be derived from considerations of the conservation of energy at a molecular level.

We can define an additional variable called the specific heat ratio, which is given the Greek symbol "gamma", which is equal to cp divided by cv:

gamma = cp / cv

Enthalpy:

Enthalpy, a property of a thermodynamic system, is the sum of the system's internal energy and the product of its pressure and volume. It is a state function used in many measurements in chemical, biological, and physical systems at a constant pressure, which is conveniently provided by the large ambient atmosphere.

If we are dealing with a gas, it is most convenient to use forms of the thermodynamics equations based on the enthalpy of the gas. From the definition of enthalpy:

h = e + p * v

Entropy:

Entropy is a scientific concept as well as a measurable physical property that is most commonly associated with a state of disorder, randomness, or uncertainty. The term and the concept are used in diverse fields, from classical thermodynamics, where it was first recognized, to the microscopic description of nature in statistical physics, and to the principles of information theory.

What is File parsing ?

File parsing means reading a file from the outer folder inside the matlab or extracting the value of a folder provided extracting the value inside the folder. The name Parse refers in matlab reading a function, script or internal code inside the memory. Some of the common commands used are (fopen, fgetl, fgets etc).

We are using this technique to read the Thermodynamic data by using file parsing to get the value of Cp, H, and S.

Specific Heat function:

  • Here I have used the same formula as mention above. To calculate the cp, h and s value.
  • For loop is for repeating these condition for all three at once cp, h and s.
  • If loop is used here because there are condition applied for finding value all 7 co-efficient of higher temperature. And value of 7 co-efficient of lower temperature.

Code:

function [cp, h, s] = coeff(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,T)

% Given conditions
R = 8.314;
for j = 1:length(T)
    if T(j) < 300
        
        cp = (a1+a2.*T+a3.*T.^2+a4.*T.^3+a5.*T.^4).*R;
        h = (a1+(a2.*T)./2+(a3.*T.^2)./3+(a4.*T.^3)./4+(a5.*T.^4)./5+(a6./T)).*R.*T;
        s = (a1.*log(T)+a2.*T+(a3.*T.^2)./2+(a4.*T.^3)./3+(a5.*T.^4)./4+a7).*R;
    else
        cp = (a8+a9.*T+a10.*T.^2+a11.*T.^3+a12.*T.^4).*R;
        h = (a8+(a9.*T)./2+(a10.*T.^2)./3+(a11.*T.^3)./4+(a12.*T.^4)./5+(a13./T)).*R.*T;
        s = (a8.*log(T)+a9.*T+(a10.*T.^2)./2+(a11.*T.^3)./3+(a12.*T.^4)./4+a14).*R;

    end
end
end

Molecular Weight function:

  • Hence we can see the thermodynamic data contains 5 species that are oxygen, carbon, nitrogen, Hydrogen, Argon.
  • These are named as elements and are stored in array.
  • Where as their atomic masses are stored in atomic mass array.
  • Here I have used For loop so that the process is to be repeated for each and every species.
  • Strcmp command is used to compare the elements and corresponding atomic mass.
  • And the molecular weight is stored in mol_wt.
  • And Hence the molecular weight is calculated.
  • Now I have to pass this function to the main code.

Code:

function [molecular_weight_cal] = mole_wt(species_name)

elements = {'H','N','O','C','A'};
atomic_mass = [1 14 16 12 40];

mol_wt = 0;
for i = 1:length(species_name)
    for j = 1:length(elements)
        if strcmp(species_name(i),elements(j))
            mol_wt = mol_wt + atomic_mass(j);
            m = j;
        end
    end

   s = str2double(species_name(i));
    if s > 1
        mol_wt = mol_wt+(atomic_mass(m)*(s-1));
    end
end
molecular_weight_cal = mol_wt
fprintf('molecular_weight_cal of %s:- ',species_name)
fprintf('%d',molecular_weight_cal)
fprintf(' ')
fprintf('amu')
fprintf('\n')
end

Main Code :

clear all
close all
clc

f1 = fopen('THERMO.dat','r');
fgetl(f1);

G_temp = str2num(fgetl(f1));

%Global temparture 

G_low = G_temp(1);
G_mid = G_temp(2);
G_high= G_temp(3);

for i = 1:3
    fgetl(f1);
end



R = 8.3145; %J/Mol*Km (units)

for j = 1:53
    
line = fgetl(f1);

A = strsplit(line,' ');
 
% species 
species =(A{1});


% Local_temparture
Low = str2num(A{end-3});
High = str2num(A{end-2});
Mid = str2num(A{end-1});

%co-efficients 
a_line = fgetl(f1);
a = findstr(a_line,'E'); %#ok<*FSTR>

A1 = a_line(1:a(1)+3);
A1 = str2num(A1); %#ok<*ST2NM>
A2 = a_line(a(1)+4:a(2)+3);
A2 = str2num(A2);
A3 = a_line(a(2)+4:a(3)+3);
A3 = str2num(A3);
A4 = a_line(a(3)+4:a(4)+3);
A4 = str2num(A4);
A5 = a_line(a(4)+4:a(5)+3);
A5 = str2num(A5);


b_line = fgetl(f1);
b = findstr(b_line,'E');

A6 = b_line(1:b(1)+3);
A6 = str2num(A6);
A7 = b_line(b(1)+4:b(2)+3);
A7 = str2num(A7);
A8 = b_line(b(2)+4:b(3)+3);
A8 = str2num(A8);
A9 = b_line(b(3)+4:b(4)+3);
A9 = str2num(A9);
A10 = b_line(b(4)+4:b(5)+3);
A10 = str2num(A10);

c_line = fgetl(f1);
c = findstr(c_line,'E');

A11 = c_line(1:c(1)+3);
A11 = str2num(A11);
A12 = c_line(c(1)+4:c(2)+3);
A12 = str2num(A12);
A13 = c_line(c(2)+4:c(3)+3);
A13 = str2num(A13);
A14 = c_line(c(3)+4:c(4)+3);
A14 = str2num(A14);

%Specific heat, Enthalpy, Entropy @ High temparture
T_H = linspace(Low,Mid,100);
  
cp = (A1+(A2.*T_H)+(A3.*T_H.^2)+(A4.*T_H.^3)+(A5.*T_H.^4))*R;
H = (A1+(A2.*T_H/2)+(A3.*(T_H.^2)/3)+(A4.*(T_H.^3)/4)+(A5.*(T_H.^4)/5)+A6.*T_H)*R;
S = (A1*log(T_H)+(A2.*T_H)+(A3.*(T_H.^2)/2)+(A4.*(T_H.^3)/3)+(A5.*(T_H.^4)/4)+A7)*R;
    
%Specific heat, Enthalpy, Entropy @ Low temparture
T_L = linspace(Mid,High,100);

cp1 = R*(A8+(A9.*T_L)+(A10.*T_L.^2)+(A11.*T_L.^3)+(A12.*T_L.^4));
H1 = R*(A8+(A9.*T_L/2)+(A10.*(T_L.^2)/3)+(A11.*(T_L.^3)/4)+(A12.*(T_L.^4)/5)+A13.*T_L);
S1 = R*(A8*log(T_L)+(A9.*T_L)+(A10.*(T_L.^2)/2)+(A11.*(T_L.^3)/3)+(A12.*(T_L.^4)/4)+A14);

% To create a folder
str = sprintf('D:\\virat\\Nasa thermodynamics data\\%s',species);
mkdir(str);
cd(str);

figure()
plot(T_H,cp)
hold on 
plot(T_L,cp1)
hold off
xlabel('Temparature')
ylabel('Specific heat') 
legend('High temparture','low Temparature')
saveas (gcf,'Specific heat')

figure()
plot(T_H,H)
hold on
plot(T_L,H1)
hold off
xlabel('Temparature')
ylabel('Enthalpy') 
legend('High temparture','low Temparature')
saveas(gcf,'Enthalpy')

figure()
plot(T_H,S)
hold on 
plot(T_L,S1)
hold off
xlabel('Temparature')
ylabel('Entropy')
legend('High temparture','low Temparature')
saveas(gcf,'Entropy')



cd ..

%Calculating Molecular weight
elements = ['H','C','N','O','AR','(S)'];
Atomic_weight = [1,12,14,16,40,0,0,0,0];

molecular_wt = 0;
for k = 1:length(species)
  [num,status] = str2num(species(k));
  if status == 0
    ind = find(elements == species(k));
    molecular_wt(k) = Atomic_weight(ind);
  else 
  molecular_wt(k-1) = num*molecular_wt(k-1);
  end
end
t_M=sum(molecular_wt);
K = sprintf('Molecular weight of %s is %f',species,t_M);
disp(K)
end

CODE EXPLANATION :

  • Starting with clear all, close all and clc
  • By using fopen command we are opening the file of THERMO.dat and 'r' repersents to read the file

fgetl: is used to read the line in the file of which we have given and for every fgetl it will be reading the next line in the file So the fgetl is used to read the line in the file

  • By using fgetl we will be reading the first line of the file which is THERMO and for next fgetl we ill be getting the second line which is the Global_temperature and we are assinging to the G_temp
  • G_low = G_temp(1)
    G_mid = G_temp(2)
    G_high= G_temp(3)
  • And we are skipping the next three lines by creating the for loop because those are the comments

  • And then in the 4TH line we have species name in the first position and then the local_temperature and we are assinging this line to the line
  • By using the strsplit we are splitting the line with the reference of the 'space'

strsplit : Split string or character vector at specified delimiter

Delimiting characters, specified as a character vector, a 1-by-n cell array of character vectors, or a 1-by-n string array. Text specified in delimiter does not appear in the output C.

Specify multiple delimiters in a cell array or a string array. The strsplit function splits str on the elements of delimiter. The order in which delimiters appear in delimiter does not matter unless multiple delimiters begin a match at the same character in str. In that case strsplit splits on the first matching delimiter in delimiter.

  • Strsplit is used to split the line and it is assigned to the variable A so the A(1) is the species name and for the local_temp we ill be taking from the end
  • Low = str2num(A{end-3})
    High = str2num(A{end-2})
    Mid = str2num(A{end-1})

str2num: Convert character array or string to numeric array

  • And then we have to read the co-efficients which are from the next line and then by using fgetl read the next line and assigned to a variable
  • To assign each value to a variable we used findstr

findstr: Find strings within other strings

  • We ill be finding the E in the string (line) and by refering that we ill be taking the separate value and assign to a variable as shown

  • a_line = fgetl(f1);
    a = findstr(a_line,'E')

    A1 = a_line(1:a(1)+3)

  • Here a will be getting the position of the E in the line and by using the position and number of numbers to left and right we ill be getting the co_efficiants and for the second number we have to start from the first number position as shown

A2 = a_line(a(1)+4:a(2)+3)

  • lly to all co_efficiants as shown
a_line = fgetl(f1);
a = findstr(a_line,'E'); %#ok<*FSTR>

A1 = a_line(1:a(1)+3);
A1 = str2num(A1); %#ok<*ST2NM>
A2 = a_line(a(1)+4:a(2)+3);
A2 = str2num(A2);
A3 = a_line(a(2)+4:a(3)+3);
A3 = str2num(A3);
A4 = a_line(a(3)+4:a(4)+3);
A4 = str2num(A4);
A5 = a_line(a(4)+4:a(5)+3);
A5 = str2num(A5);


b_line = fgetl(f1);
b = findstr(b_line,'E');

A6 = b_line(1:b(1)+3);
A6 = str2num(A6);
A7 = b_line(b(1)+4:b(2)+3);
A7 = str2num(A7);
A8 = b_line(b(2)+4:b(3)+3);
A8 = str2num(A8);
A9 = b_line(b(3)+4:b(4)+3);
A9 = str2num(A9);
A10 = b_line(b(4)+4:b(5)+3);
A10 = str2num(A10);

c_line = fgetl(f1);
c = findstr(c_line,'E');

A11 = c_line(1:c(1)+3);
A11 = str2num(A11);
A12 = c_line(c(1)+4:c(2)+3);
A12 = str2num(A12);
A13 = c_line(c(2)+4:c(3)+3);
A13 = str2num(A13);
A14 = c_line(c(3)+4:c(4)+3);
A14 = str2num(A14);
  • We have to write a formula for low and high temperatures The FIRST 7 coefficients are HIGH-temperature coefficients and the SECOND 7 coefficients are LOW-temperature coefficients.
%Specific heat, Enthalpy, Entropy @ High temparture
T_H = linspace(Low,mid,100);
  
cp = (A1+(A2.*T_H)+(A3.*T_H.^2)+(A4.*T_H.^3)+(A5.*T_H.^4))*R;
H = (A1+(A2.*T_H/2)+(A3.*(T_H.^2)/3)+(A4.*(T_H.^3)/4)+(A5.*(T_H.^4)/5)+A6.*T_H)*R;
S = (A1*log(T_H)+(A2.*T_H)+(A3.*(T_H.^2)/2)+(A4.*(T_H.^3)/3)+(A5.*(T_H.^4)/4)+A7)*R;
    
%Specific heat, Enthalpy, Entropy @ Low temparture
T_L = linspace(mid,High,100);

cp1 = R*(A8+(A9.*T_L)+(A10.*T_L.^2)+(A11.*T_L.^3)+(A12.*T_L.^4));
H1 = R*(A8+(A9.*T_L/2)+(A10.*(T_L.^2)/3)+(A11.*(T_L.^3)/4)+(A12.*(T_L.^4)/5)+A13.*T_L);
S1 = R*(A8*log(T_L)+(A9.*T_L)+(A10.*(T_L.^2)/2)+(A11.*(T_L.^3)/3)+(A12.*(T_L.^4)/4)+A14);
  • And then we have to create a folder to save the plots in thier respective folders to do that we will be using the mkdir

mkdir: To create a folder

cd: Change current folder

  • str = sprintf('D:viratNasa thermodynamics data%s',Spices);
    mkdir(str);
    cd(str);
  • plotting the enthalpy, entropy and specific heats for the low and high temperatures co_efficiants
figure()
plot(T_H,cp)
hold on 
plot(T_L,cp1)
hold off
xlabel('Temparature')
ylabel('Specific heat') 
legend('High temparture','low Temparature')
saveas (gcf,'Specific heat')

figure()
plot(T_H,H)
hold on
plot(T_L,H1)
hold off
xlabel('Temparature')
ylabel('Enthalpy') 
legend('High temparture','low Temparature')
saveas(gcf,'Enthalpy')

figure()
plot(T_H,S)
hold on 
plot(T_L,S1)
hold off
xlabel('Temparature')
ylabel('Entropy')
legend('High temparture','low Temparature')
saveas(gcf,'Entropy')
  • after every plot saveas will save the plot as fig in the respective folder

saveas: Save figure to specific file format gcf is the handle to save the plot as image

  • And the we have to give the cd.. command to come out from the particular folder

cd.. : Change current folder

../ One level up from the current folder.
Multiples of ../ Multiple levels up from the current folder.
./ A path relative to the current folder. cd assumes that the path is relative to the current folder even if the ./ is not present.
  • Calculating the molecular weight
%Calculating Molecular weight
elements = ['H','C','N','O','AR','(S)'];
Atomic_weight = [1,12,14,16,40,0,0,0,0];


for k = 1:length(species)
  [num,status] = str2num(species(k));
  if status == 0
    ind = find(elements == species(k));
    molecular_wt(k) = Atomic_weight(ind);
  else 
  molecular_wt(k-1) = num*molecular_wt(k-1);
  end
end
t_M=sum(molecular_wt)
K = sprintf('Molecular weight of %s is %f',species,t_M);
disp(K)
  • Create two strings of atom symbol and the atomic weight as shown
  • elements = ['H','C','N','O','AR','(S)'];
    Atomic_weight = [1,12,14,16,40,0,0,0,0];
  • Create a for loop to the length of the species convert the species to number by using str2num command with 2 results
  • In species if it is number status will be 1 it is true if it is an letter(element) the status will be 0
  • By applying if condition status == 0
  • ind = find(elements == species(k)) we are finding the elements in species in the elements array by using find command
  • Then the molecular_wt = atomic_weight(ind) else
  • molecular_wt(k-1) is it will take the previous position of the species equals to number in the species multiply to the molecular_wt(k-1) so that the previous position to the number in the species will be changed
  • And then display the molecular_wt

RESULTS

  • molecular weight of each species display it in the command window.

  • Save the plots as images with appropriate names and dump them in separate folders for each species with suitable name is done AUTOMATICALLY
  • Plots for o2, N2 and Co2 species

PLOTS OF O2

specific heats(cp) enthalpy(H) entropy(S)

PLOTS OF N2

specific heats(cp) enthalpy(H) entropy(S)

PLOTS OF CO2

specific heats(cp) enthalpy(H) entropy(S)

ERRORS FACED

  • Forgot to give double slash

  • After CH2 it is CH2(S) so it shows that there are differant number of elements

Conclusion:

  • The thermodynamic data has been parsed, and thermodynamic data has been calculated and plotted for each species.
  • It has been clearly seen that the data of cp, h and s vary at different temperature range from low-temperature to high-temperature at each compound.
  • Hence we have calculated the Molecular weight of different species and stored the data using file parsing technique.

References:

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.