Challenge

Project 1 - Parsing NASA thermodynamic data

SShubham Ajitsaria30 Nov 2020

AIM: To write a function that extracts the 14 coefficients and calculates the enthalpy, entropy, and specific heats for all the species in the data file. Also, Calculate the molecular weight of each species and display it in the command window.

Theory:

a: Parsing

Parsing a file means reading in a data stream of some sort and building an in-memory model of the semantic content of that data. This aims at facilitating performing some kind of transformation on the data.

b: Specific heat

In Thermodynamics it is the ratio of the quantity of heat required to raise the temperature of a body one degree to that required to raise the temperature of an equal mass of water one degree.

Where Cp=Specific Heat

R= Universal Gas Constant

C: Enthalpy

Enthalpy (H) is the sum of the internal energy (U) and the product of pressure and volume (PV) given by the equation: H=U+PV. When a process occurs at constant pressure, the heat evolved (either released or absorbed) is equal to the change in enthalpy.

d: Entropy

It is the level of randomness (or disorder) of a system. It could also be thought of as a measure of the energy dispersal of the molecules in the system. Microstates are the number of different possible arrangements of molecular position and kinetic energy at a particular thermodynamic state.

Main Code:

clear all
close all
clc
%Opening and reading the file
f = fopen('THERMO.dat','r')

%Header
line1 = fgetl(f)

%Global Temperature data
line2 = fgetl(f)

%1st 3 lines are comments
com1 = fgetl(f)
com2 = fgetl(f)
com3 = fgetl(f)
A = strsplit(line2,' ');
global_low_temp = str2num(A{2})
global_mid_temp = str2num(A{3})
global_high_temp = str2num(A{4})

%Species/Elements information
for  i =1:53
    species = fgetl(f);
    species_name = species(1:7)
    t1 = fgetl(f);
    l1 = strfind(t1,'E');
    c1 = str2num(t1(1:(l1(1)+3)));
    c2 = str2num(t1((l1(1)+4):(l1(2)+3)));
    c3 = str2num(t1((l1(2)+4):(l1(3)+3)));
    c4 = str2num(t1((l1(3)+4):(l1(4)+3)));
    c5 = str2num(t1((l1(4)+4):(l1(5)+3)));
    
    t2 = fgetl(f);
    l2 = strfind(t2,'E')
    
    c6 = str2num(t2(1:(l2(1)+3)));
    c7 = str2num(t2((l2(1)+4):(l2(2)+3)));
    c8 = str2num(t2((l2(2)+4):(l2(3)+3)));
    c9 = str2num(t2((l2(3)+4):(l2(4)+3)));
    c10 = str2num(t2((l2(4)+4):(l2(5)+3)));
    
    t3 = fgetl(f);
    l3 = strfind(t3,'E');
    
    c11 = str2num(t3(1:(l3(1)+3)));
    c12 = str2num(t3((l3(1)+4):(l3(2)+3)));
    c13 = str2num(t3((l3(2)+4):(l3(3)+3)));
    c14 = str2num(t3((l3(3)+4):(l3(4)+3)));
    
    %Creating folder for each species
    mkdir(['D:\USERS D\SHUBHAM D\DOWNLOADS D\designs save\SKILL LYNC PROJECTS\Matlab\Test program\nasa data',species_name]);
    
    
    
    %Temperature ranges assumed for calculation perpose
    temp = linspace(300,5000,100);
    
    %Universal gas constant
    R = 8.314;
    
    %calculation of properties based on lower and upper temperature
    for j = 1:length(temp)
        if temp(j)<3000
            Cp = (c1+c2.*temp+c3.*temp.^2+c4.*temp.^3+c5.*temp.^4).*R;
            H = (c1+(c2.*temp)./2+(c3.*temp.^2)./3+(c4.*temp.^3)./4+(c5.*temp.^4)./5+(c6./temp)).*R.*temp;
            S = (c1.*log(temp)+c2.*temp+(c3.*temp.^2)./2+(c4.*temp.^3)./3+(c5.*temp.^4)./4+c7).*R;
        else
            Cp = (c8+c9.*temp+c10.*temp.^2+c11.*temp.^3+c12.*temp.^4).*R;
            H = (c8+(c9.*temp)./2+(c10.*temp.^2)./3+(c11.*temp.^3)./4+(c12.*temp.^4)./5+(c13./temp)).*R.*temp;
            S = (c8.*log(temp)+c9.*temp+(c10.*temp.^2)./2+(c11.*temp.^3)./3+(c12.*temp.^4)./4+c14).*R;
        end
    end
    
    %plotting
    
    %temperature vs Entropy
    figure(1)
    plot(temp,S,'linewidth',4,'color','b')
    title(species_name)
    xlabel('Temperature')
    ylabel('Entropy')
    
    %temperature vs Enthalpy
    figure(2)
    plot(temp,H,'linewidth',4,'color','k')
    title(species_name)
    xlabel('Temperature')
    ylabel('Enthalpy')
    
    %temperature vs Specific heat
    figure(3)
    plot(temp,Cp,'linewidth',4,'color','y')
    title(species_name)
    xlabel('Temperature')
    ylabel('Specific heat')
    directory =(['D:\USERS D\SHUBHAM D\DOWNLOADS D\designs save\SKILL LYNC PROJECTS\Matlab\Test program\nasa data',species_name]);
    %changing the directory
    cd(directory)
    
 
   %saving the plots
    
    saveas(1,'Entropy.jpg');
    saveas(2,'Enthalpy.jpg');
    saveas(3,'Specific heat.jpg');
    %switching back to orignal directory
    cd('D:\USERS D\SHUBHAM D\DOWNLOADS D\designs save\SKILL LYNC PROJECTS\Matlab\Test program\nasa data')
    M = molecular_mass(species_name)
    
end
fclose(f)

Code Flow Steps Explanations

Used fopen command to open the THERMO.DAT file to read

Used fgetl to get header comments and global temperature data and stored in different variables.

strsplit to split the provide global temperature which was in the form of combined data

str2num to convert string to numbers. to extract temperature and store in other variables

Using for loop to extract information of each species by repeating the extraction algorithm.

strfind to find separater between constants and combine exponential part

mkdir was used to create a folder directory to store the result data

if-else is used to solve for low temp and high temp separately based on the value of temp.

Also common commands such as plot (to plot graph) and cd (to change directory was used)

molecular_mass function is formed so as to calculate the molecular weight of the given elements

Molecular weight calculation code:

function M = molecular_mass(species_name)
species = upper(species_name)
constituents = ['H','C','N','O','A'];
atomic_weight = [1.00749 12.0107 14.0067 15.9994 39.948];
M = 0;
for i=1:length(species)
    for j=1:length(constituents)
        if strcmp(species(i),constituents(j))
            M = M+atomic_weight(j);
            u=j;
        end
    end
    n =str2num(species(i));
    if n>1
     M = M + atomic_weight(u)*(n-1);
    end
end
end

Result:

The plot of Element O2 (in order of Enthalpy, Entropy, Specific heat)