NASA File parsing to calculate thermodynamics properties using MATLAB

---

NASA came up with polynomials that can be used 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:

• Write the function that extracts the 14 coefficients and calculates the enthaply, entropy and specific heat for all the species in the data file
• Calculate the molecular weight of each species and display it in the command window
• Plot the Cp, Enthaply and entropy for local temperature range(low temperature :high temperature) specific for each species
• Save the plots as images with appropriate names and dump them in separate folder for each species with suitable name. All these should be generated automatically

Procedure:

In this project, We will be reading NASA thermodynamics data and the calculate the theromodynamics properties such as Specific heat, enthaply, entropy and molecular weight of various species and plot the respective graph. The required data is stored in the file named \'THERMO.dat\' and data is used for the project. The formula used for calculating the value Cp, Enthaply and entropy are given below:

`(Cp)/R= a1 + a2 T+a3T^2 +a4T^3+a5T^4`

`H/(RT)=a1 +a2T/2+a3T^2/3+a4T^3/4+a5T^4/5+(a6)/T`

`S/R=a1lnT+a2T+a3T^2/2+a4T^3/3+a5T^4/4+a5T^4/4+a7`

Here.

Cp - Specific heat

H- Enthaply

S- Entropy

T- Temperature

a1,a2,a3,a4,a5,a6,a7 are the coefficients

The first 7 coefficients are low temperature coefficients and the second 7 coefficients are high temperature coefficients. The mid temparture range is 100K and if anything goes higher than that, its considered as high temperature.

STEPS FOR CODE

• Intially, Open the file THERMO.dat while using fopen command
• fgetl command is used to to read the first line of the data file, after fscan is used used for to read global low, mid and high temp values are used
• strfind command is used to find the position of \'E\' and in line to seperate coefficient
• New directory is made usinng mkdir command to create a folder with the species name
• To calculate Entropy, Enthalpy and specific Heat functions are defined
• after then graph has been plotted T vs Cp, T vs H, T vs S
• using cd directory is changed and the molecular mass is calculated using function which has been shown below

clear all 
close all
clc
%universal gas constant
R=8.314

f1=fopen(\'THERMO.dat\',\'r\')           %open thermo.dat file from the same directory where codes has been saved
f2=fopen(\'Molecular mass.txt\',\'w\')   %open notepad save molecular mass for all species and write inside the notepad
line1=fgetl(f1)                      % fgetl command is used to read the fiorst line of the file
%intial lines
for i=1:4
    fgetl(f1)
end
  
T=linspace(300,1000,5000)              

for i=1:52
    line=fgetl(f1)                    %extracting first and species name
    A=strsplit(line,\' \')
    species=A(1)
    
    line2=fgetl(f1)
    sf=strfind(line2,\'E\')            %strfind- to find string
    
    %extracting higher temp coeficients
    cf1=line2(1:sf(1)+3)
    cf2=line2(sf(1)+4:sf(2)+3)
    cf3=line2(sf(2)+4:sf(3)+3)
    cf4=line2(sf(3)+4:sf(4)+3)
    cf5=line2(sf(4)+4:sf(5)+3)
    
    line3=fgetl(f1)
    cf6=line3(1:sf(1)+3)
    cf7=line3(sf(1)+4:sf(2)+3)
    
    %extracting lower temp coeeficients
    cfe1=line3(sf(2)+4:sf(3)+3)
    cfe2=line3(sf(3)+4:sf(4)+3)
    cfe3=line3(sf(4)+4:sf(5)+3)
    
    line4=fgetl(f1)
    
    cfe4=line4(1:sf(1)+3)
    cfe5=line4(sf(1)+4:sf(2)+3)
    cfe6=line4(sf(2)+4:sf(3)+3)
    cfe7=line4(sf(3)+4:sf(4)+3)
   
    %temperature coefficient
    hc1=str2num(cf1);
    hc2=str2num(cf2);
    hc3=str2num(cf3);
    hc4=str2num(cf4);
    hc5=str2num(cf5);
    hc6=str2num(cf6);
    hc7=str2num(cf7);
    
    lc1=str2num(cfe1);
    lc2=str2num(cfe2);
    lc3=str2num(cfe3);
    lc4=str2num(cfe4);
    lc5=str2num(cfe5);
    lc6=str2num(cfe6);
    lc7=str2num(cfe7);
    %low temp cofficient
    lc=[lc1 lc2 lc3 lc4 lc5 lc6 lc7];
    %high temp coefficient
    hc=[hc1 hc2 hc3 hc4 hc5 hc6 hc7];
    
    for j=1:length(T)
        t=1000;
        if T(j)<=t
            Cp(j)=R*(lc(1)+lc(2)*T(j)+lc(3)*(T(j))^2+lc(4)*(T(j))^3+lc(5)*(T(j))^4);
            H(j)=R*T(j)*(lc(1)+lc(2)*(T(j)/2)+lc(3)*(((T(j))^2)/3)+lc(4)*(((T(j))^3)/4)+lc(5)*(((T(j))^4)/5)+lc(6)/T(j));
            S(j)=R*(lc(1)*log(T(j))+lc(2)*T(j)+lc(3)*((T(j))^2/2)+lc(4)*((T(j))^3/3)+lc(5)*((T(j))^4/4)+lc(7));
        else
            Cp(j)=R*(hc(1)+hc(2)*T(j)+hc(3)*(T(j))^2+hc(4)*(T(j))^3+hc(5)*(T(j))^4);
            H(j)=R*T(j)*(hc(1)+hc(2)*(T(j)/2)+hc(3)*(((T(j))^2)/3)+hc(4)*(((T(j))^3)/4)+hc(5)*(((T(j))^4)/5)+hc(6)/T(j));
            S(j)=R*(hc(1)*log(T(j))+hc(2)*T(j)+hc(3)*((T(j))^2/2)+hc(4)*((T(j))^3/3)+hc(5)*((T(j))^4/4)+hc(7));
        end
    end
    
    spec=char(species)
    fd=sprintf(\'%s\',spec);
    cur_dir=pwd;
    mkdir(spec);
    cd(spec);
    
    %plot the first graph Temp Vs Soecific heat
    figure(1)
    plot(T,Cp);
    xlabel(\'Temperature (K)\');
    ylabel(\'Specific heat (KJ)\');
    title(\'Specific Heat Vs Temperature\');
    
    %plot the second graph temp Vs Enthaply
    figure(2)
    plot(T,H)
    xlabel(\'Temperature (K)\');
    ylabel(\'Enthalpy (KJ/mol)\');
    title(\'Enthalpy Vs Temperature\');
    
    %Plot the third graph Temperature Vs Entropy
    figure(3)
    plot(T,S)
    xlabel(\'Temperature (K)\');
    ylabel(\'Entropy (KJ/mol-K)\');
    title(\'Entropy Vs Temperature\');
    
    %saving the plot that got generated
    saveas(figure(1),\'Specific_heat.png\')
    saveas(figure(2),\'Enthaply.png\')
    saveas(figure(3),\'Entropy.png\')
    
    cd(cur_dir);
    
    %Molecular weight of particular species
    mm(i)=Molecular_weights(spec);
    fprintf(f2,\'%s \\n %s %f \\n\',spec,\'Molecular mass=\',mm(i))
end
fclose(f2);  %close file F2
fclose(f1)    %close file f1

    

Function for Molecular weights:

%Function for calculating molecular weight
function M=Molecular_weights(species_element)
atoms=[\'H\',\'C\',\'N\',\'O\',\'R\'];
disp(\'Molecular weight of \')
disp(species_element);
disp(\':\');
Atomic_weight=[1 12 14 16 40]
M=0;
for i=1:length(species_element)
    for j=1:length(atoms)
        if strcmp(species_element(i),atoms(j))
            M=M+Atomic_weight(j)
            p=j;
        end
    end
    n=str2num(species_element(i));
    if n>1
        M=M+(Atomic_weight(p))*(n-1);
    end
end
fprintf(\'%d\',M)
disp(\'\');
end

For N2- Plots are shown below

For O2 Plots are shown below

For Co2 Plots are shown below

Folder Screenshot where all files saved:

Molecular mass for all species Screenshot:

Sharable link for the molecular wewights for all the species:

https://drive.google.com/open?id=1mTf6kkaHXATF4GzWm_FzzCkD_HWdscc2

Reference:

https://in.mathworks.com/matlabcentral/answers/196582-best-way-to-parse-text-file

https://in.mathworks.com/matlabcentral/answers/47477-parsing-formatted-text-file-quickly