Challenge

Project 1 - Parsing NASA thermodynamic data

RRhitwik Sundar9 Mar 2021

Aim:

Parsing NASA thermodynamic data to calculates the enthalpy, entropy and specific heats for all the species in the given data file.

Objective:

Write a code to calculate the Cp,S and H of all the species in the in the given data file by using the 14 co-efficents, like for first 7 coefficients are high-temperature coefficients and the second 7 coefficients are low-temperature coefficients.Use local temperature for every species.
Calculate the molecular weight of each species and display it in the command window.
Make a code to save the plots of every species in different folders automatically.

Abstract:

In this project we are coding for NASA thermodynamic data about different species.In each species having 14 coefficients and these coefficients are given to the formula to calculate Specific heat, Enthalpy, Entropy. In each of these formula having coefficients according to temperature range coefficient need to assign on each formulas and also for each species Enthalpy ,Entropy and Specific heat will be calculating.

Enthalpy ,Entropy and Specific heats for every species need to plotted and this plots should saved in a folder with species name and final we need to calculate the molecular weight of each species and display in the results.

Governing Equation:

Specific Heat(Cp):Based on thermodynamics ,the specific heat capacity of a substance is the heat capacity of a sample of the substance divided by the mass of the sample. In layman's words,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 temperature. The SI unit of specific heat is joule per kelvin and kilogram J/(K kg).

The above formula is which we use to calculate Specific heat(Cp) in this project.

Where,R is the universal gas constant

T is the temperature

a1 ,a2 ,a3 ,a4 ,a5 these are the coefficients of the temperature

Enthalpy: It is concept used in science and engineering when the heat and work need to be calculated. if a substance changes at constant pressure at this juncture enthalpy tells how much heat and work was added or removed from the substance . Enthalpy is analogous to energy but not the similar. During the process of substance enlarging and shrinking ,energy is consumed or released.

The given below formula that we use to calculate enthalpy in this project,

Where,R is the universal gas constant

T is the temperature

a1 ,a2 ,a3 ,a4 ,a5 , a6 these are the coefficients of the temperature.

Entropy:defined as measure of the amount of energy which is unavailable to do work .it is a measure of the number of possible arrangements of the atom in a system that can have ,so it could be called as measure of randomness.

The given below formula that we use to calculate Entropy in this project,

Where,R is the universal gas constant

T is the temperature

a1 ,a2 ,a3 ,a4 ,a5 , a7 these are the coefficients of the temperature.

Code for parsing NASA thermodynamic data:

clear all
close all
clc

%To open the data file 
data = fopen('THERMO.dat','r');     %opening the data file for reading the data
fgetl(data);
fgetl(data);
fgetl(data);
fgetl(data);
fgetl(data);

for j=1:53        %for 53 species
    

   
    %Header information
    A =fgetl(data);

    %splitting the line

    A =strsplit(A,' ');               %For splitting the next line of the data with space in between

    %Reading Global Temperatures

  Low_temperature =str2num(A{end-3});    %To display low temperature from the data
  Medium_temperature =str2num(A{end-1}); %To display medium temperature from the data
  High_temperature =str2num(A{end-2});   %To display high temperature from the data

    %To display species name from the data

     species_name= (A{1});

    %To read the  1set of coefficients
    L1= fgetl(data);
    a = findstr(L1,'E');
    a1= L1(1:a(1)+3);
    a1= str2num(a1);
    a2= L1(a(1)+4:a(2)+3);
    a2= str2num(a2);
    a3= L1(a(2)+4:a(3)+3);
    a3= str2num(a3);
    a4= L1(a(3)+4:a(4)+3);
    a4= str2num(a4);
    a5= L1(a(4)+4:a(5)+3);
    a5= str2num(a5);

    %To read the 2nd set of coefficients
    L2= fgetl(data);
    b = findstr(L2,'E');
    a6= L2(1:b(1)+3);
    a6= str2num(a6);
    a7= L2(b(1)+4:b(2)+3);
    a7= str2num(a7);
    a8 = L2(b(2)+4:b(3)+3);
    a8 = str2num(a8);
    a9 = L2(b(3)+4:b(4)+3);
    a9 = str2num(a9);
    a10 = L2(b(4)+4:b(5)+3);
    a10 = str2num(a10);

    %To read the 3rd set of coefficients
    L3 =fgetl(data);
    c  = findstr(L3,'E');
    a11= L3(1:c(1)+3);
    a11= str2num(a11);
    a12= L3(c(1)+4:c(2)+3);
    a12=str2num(a12);
    a13=L3(c(2)+4:c(3)+3);
    a13=str2num(a13);
    a14=L3(c(3)+4:c(4)+3);
    a14=str2num(a14);
    
    %Reading the temperature array
     Temperature=linspace(Low_temperature,High_temperature,1000);  %Creating a array  of 1000 temperature values between low temperature to high temperature
i = 1;
     %Universal gas constant value
     R= 8.31445;
        for  T=linspace(Low_temperature,High_temperature,1000)
          if (T>Medium_temperature)
              %calculate Cp ,H ,S
          
             Cp(i)=((a1)+(a2*T)+(a3*T^2)+(a4*T^3)+(a5*T^4))*R;

             H(i)=((a1)+(a2*T/2)+(a3*T^2/3)+(a4*T^3/4)+(a5*T^4/5)+(a6/T))*R*T;

             S(i)=((a1*log(T))+(a2*T)+(a3*T^2/2)+(a4*T^3/3)+(a5*T^4/4)+(a7))*R;
          else
              
              % Formula to calculate Cp ,H ,S
    
             Cp(i)=((a8)+(a9*T)+(a10*T^2)+(a11*T^3)+(a12*T^4))*R;

             H(i)=((a8)+(a9*T/2)+(a10*T^2/3)+(a11*T^3/4)+(a12*T^4/5)+(a13/T))*R*T;

             S(i)=((a8*log(T))+(a9*T)+(a10*T^2/2)+(a11*T^3/3)+(a12*T^4/4)+(a14))*R;   
          end
          %incrementing each values
i = i+1;
          
         
        end
    
   %Creating folder for each species      
    mkdir (['C:\Users\Skill_Lync\Desktop\matlab\NASA_thermodynamic_data\Species_',species_name]);
    
    %Plotting the results
    
    figure(1)
    plot(Temperature,Cp,'linewidth',2,'color','r')
    title('Specific Heat vs Temperature')
    xlabel('Temperature')
    ylabel('Specific Heat')
    savings=(['C:\Users\Skill_Lync\Desktop\matlab\NASA_thermodynamic_data\Species_',species_name]);
    cd(savings);
    saveas(1,'Specific_heat')
   
    figure(2)
    plot(Temperature,H,'linewidth',2,'color','g')
    title('Enthalpy vs Temperature')
    xlabel('Temperature')
    ylabel('Enthalpy')
      savings=(['C:\Users\Skill_Lync\Desktop\matlab\NASA_thermodynamic_data\Species_',species_name]);
    cd(savings);
    saveas(2,'Enthalpy')
    
    figure(3)
    plot(Temperature,S,'linewidth',2,'color','b')
    title('Entropy vs Temperature')
    xlabel('Temperature')
    ylabel('Entropy')
      savings=(['C:\Users\Skill_Lync\Desktop\matlab\NASA_thermodynamic_data\Species_',species_name]);
    cd(savings);
    saveas(3,'Entropy')
       
    %finding Molecular mass
    Elements=['H','C','O','N','A'];
    atomic_masses= [1.00784,12.0107,15.9994,14.0067,39.948];
    %Calculating molecular mass for different species
    Molecular_weight=0;
    for r=1:length(species_name)
        [num,status]=str2num(species_name(r));
        if (status)
            Molecular_weight= Molecular_weight+(atomic_masses(n)*(num-1));
        else
            for k=1:length(Elements)
                
                if species_name(r)==Elements(k)
                    Molecular_weight= Molecular_weight+atomic_masses(k);
                    n=k;
                end
            end
        end
        
       
    end  
    
      fprintf("Molecular_weight of %s is :%f \n",species_name,Molecular_weight);  
   


end

Step wise explanation

clear all
close all
clc

This command is used to clear all the data in workspace ,command window and figures.

data = fopen('THERMO.dat','r');

This code is used for opening the data and r is used to read the data.

fgetl(data);
fgetl(data);
fgetl(data);
fgetl(data);
fgetl(data);

This code used for to skip the first 5 line which is not neccesary for our project.

for j=1:53

This command is used for a loop from 1:53 as we have 53 species in total.

 A =fgetl(data);

This command is used to get the 1st line of data we need beacuse we had negelect the first 5 lines so inorder to get the 6th line we use this code.

 A =strsplit(A,' ');

This code is used for splitting the next line of the data with space in between,for splitting the data strsplit is using and inorder to get a space " " is used.

Low_temperature =str2num(A{end-3});   
Medium_temperature =str2num(A{end-1});   
High_temperature =str2num(A{end-2});

To read the temperatures from the data we use this command as we know the last three cells of species are the temperature , in here str2num is used to convert string to number and the values is subtracted by the end values.

     species_name= (A{1});

To display species name from the data and all the species name will save in species_name.

  L1= fgetl(data);
    a = findstr(L1,'E');
    a1= L1(1:a(1)+3);
    a1= str2num(a1);
    a2= L1(a(1)+4:a(2)+3);
    a2= str2num(a2);
    a3= L1(a(2)+4:a(3)+3);
    a3= str2num(a3);
    a4= L1(a(3)+4:a(4)+3);
    a4= str2num(a4);
    a5= L1(a(4)+4:a(5)+3);
    a5= str2num(a5);

We need to read the coefficient so it starts from the 7th line so we use 1 more fgetl is used ,in next line we using findstr to find the string 'E' each coloumn has a 'E' so we find that E from the numbers and after that E ,there is three values (example 3.28253784E+00) here are three values after the 'E') in every columns it will be like this, a(1)+3 is used to add up the next 3 values. then this value is converted into a number by using str2num code.

For 2nd coefficients a(1)+4 is used and 2nd coefficient value end to a(2)+3 from this step 2nd coefficient values will be assign in a(2).

As the same logic used to find 2nd and 3rd set of coefficients.

    %To read the 2nd set of coefficients
    L2= fgetl(data);
    b = findstr(L2,'E');
    a6= L2(1:b(1)+3);
    a6= str2num(a6);
    a7= L2(b(1)+4:b(2)+3);
    a7= str2num(a7);
    a8 = L2(b(2)+4:b(3)+3);
    a8 = str2num(a8);
    a9 = L2(b(3)+4:b(4)+3);
    a9 = str2num(a9);
    a10 = L2(b(4)+4:b(5)+3);
    a10 = str2num(a10);

So now we take the 8th line and take the next 5coefficients,the above same steps are used for reading the set of coefficients.

    %To read the 3rd set of coefficients
    L3 =fgetl(data);
    c  = findstr(L3,'E');
    a11= L3(1:c(1)+3);
    a11= str2num(a11);
    a12= L3(c(1)+4:c(2)+3);
    a12=str2num(a12);
    a13=L3(c(2)+4:c(3)+3);
    a13=str2num(a13);
    a14=L3(c(3)+4:c(4)+3);
    a14=str2num(a14);

Now we take the 9th line and take the next 5 coefficients .

Temperature=linspace(Low_temperature,High_temperature,1000);

Creating an array of 1000 temperature values between low temperature to high temperature.

R= 8.31445 it is the universal gas constant value .

 for  T=linspace(Low_temperature,High_temperature,1000)
          if (T>Medium_temperature)
              %calculate Cp ,H ,S
              
             %High temperature coefficients
             Cp(i)=((a1)+(a2*T)+(a3*T^2)+(a4*T^3)+(a5*T^4))*R;

             H(i)=((a1)+(a2*T/2)+(a3*T^2/3)+(a4*T^3/4)+(a5*T^4/5)+(a6/T))*R*T;

             S(i)=((a1*log(T))+(a2*T)+(a3*T^2/2)+(a4*T^3/3)+(a5*T^4/4)+(a7))*R;
          else
              
             % Low temperature coefficients
             Cp(i)=((a8)+(a9*T)+(a10*T^2)+(a11*T^3)+(a12*T^4))*R;

             H(i)=((a8)+(a9*T/2)+(a10*T^2/3)+(a11*T^3/4)+(a12*T^4/5)+(a13/T))*R*T;

             S(i)=((a8*log(T))+(a9*T)+(a10*T^2/2)+(a11*T^3/3)+(a12*T^4/4)+(a14))*R;   
          end

As in the question The First 7 coefficients are High-temperature coefficients and the Second 7 coefficients are LOW-temperature coefficients. so we are using the if command.

Cp(i)=((a1)+(a2*T)+(a3*T^2)+(a4*T^3)+(a5*T^4))*R;

H(i)=((a1)+(a2*T/2)+(a3*T^2/3)+(a4*T^3/4)+(a5*T^4/5)+(a6/T))*R*T;

S(i)=((a1*log(T))+(a2*T)+(a3*T^2/2)+(a4*T^3/3)+(a5*T^4/4)+(a7))*R;

This command is used to find the Cp ,H , S .

 mkdir (['C:\Users\Skill_Lync\Desktop\matlab\NASA_thermodynamic_data\Species_',species_name]);

This is the used to creating folder for each species.

    figure(1)
    plot(Temperature,Cp,'linewidth',2,'color','r')
    title('Specific Heat vs Temperature')
    xlabel('Temperature')
    ylabel('Specific Heat')
    savings=(['C:\Users\Skill_Lync\Desktop\matlab\NASA_thermodynamic_data\Species_',species_name]);
    cd(savings);
    saveas(1,'Specific_heat')

This command is used to plot the figures , so it will plot temperature vs Cp with a linewidth of 2,title will give specific heat vs temperature as title,xlabel will give Temperature as xlabel and y label will give Specific heat as ylabel.

savings(['C:\Users\Skill_Lync\Desktop\matlab\NASA_thermodynamic_data\Species_',species_name]) it will give a path to save the plots where i needed.

saveas(1,'Specific_heat) it will save by the name.

 figure(2)
    plot(Temperature,H,'linewidth',2,'color','g')
    title('Enthalpy vs Temperature')
    xlabel('Temperature')
    ylabel('Enthalpy')
      savings=(['C:\Users\Skill_Lync\Desktop\matlab\NASA_thermodynamic_data\Species_',species_name]);
    cd(savings);
    saveas(2,'Enthalpy')
    
    figure(3)
    plot(Temperature,S,'linewidth',2,'color','b')
    title('Entropy vs Temperature')
    xlabel('Temperature')
    ylabel('Entropy')
      savings=(['C:\Users\Skill_Lync\Desktop\matlab\NASA_thermodynamic_data\Species_',species_name]);
    cd(savings);
    saveas(3,'Entropy')

These code is same as figure(1) only Cp ,H, S is changing all the rest of them are same as above.

 Elements=['H','C','O','N','A'];
    atomic_masses= [1.00784,12.0107,15.9994,14.0067,39.948];
    %Calculating molecular mass for different species
    Molecular_weight=0;
    for r=1:length(species_name)
        [num,status]=str2num(species_name(r));
        if (status)
            Molecular_weight= Molecular_weight+(atomic_masses(n)*(num-1));
        else
            for k=1:length(Elements)
                
                if species_name(r)==Elements(k)
                    Molecular_weight= Molecular_weight+atomic_masses(k);
                    n=k;
                end
            end
        end
        
       
    end  
    
      fprintf("Molecular_weight of %s is :%f \n",species_name,Molecular_weight);  
   


end

This command is used to calculate molecular mass of different species, so elements and and their corresponding element atomic masses are assigned in a array.Next a for loop is defined to run all the species molecular weight,inside the for loop looking for the status from the str2num(species_name)the if or else command is taken to calculating the molecular weight ,in for loop if the species having number it will save on num. or if it is a letter it will update on status .

For an example if we are taking 'H2O' the molecular weight of H is calculated from the else command so what is happening here is the num value will be [] and status will be 0 , by taking next '2' num will be 2 and status will be 1 so the first if command will go on and taking next charcter 'O' num will be [] and status will be 0 so it will go to next else command ,So all the values of 'H2O' is update in Molecular_weight and it will print by fprintf.

So all the species calculated in the for loop .

Errors i faced

While doing this program i faced a lot of error regarding of the assigning of incorrect variables, and declaration of inputs in for loops.i had give wrong inputs while finding molecular weight of species . by looking each of errors i do it correctly and understand the concepts and keep doing the projects.

And lot of errors i do face,by keep checking where the error comes i do cleared all the errors and make the code success.

Results:

plots as images with appropriate names and save automatically them in separate folders for each species with suitable name

Plots for O2, N2 and Co2 species:

Plot for O2

Cp Vs Temperature(O2)

Enthalpy vs Temperature(O2)

Entropy vs Temperature(O2)

Plot for N2

Cp Vs Temperature(N2)

Enthalpy vs Temperature(N2)

Entropy vs Temperature(N2)

Plot for Co2

Cp Vs Temperature(CO2)

Enthalpy vs Temperature(CO2)

Entropy vs Temperature(CO2)

Molecular weight of each species:

Conclusion

Hence we had done the Parsing of the NASA file with plot generation,plot for 53 species are generated and save it automatically to a folder by their species name and also 3 plots (Cp, H, S) were included in these folders for different species.

And finally we had find 53 species molecular weight and generated in the command window ,hence we had done all the required details for this project.

Project 1 - Parsing NASA thermodynamic data

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.