NASA Thermodynamic parsing by MAT LAB

AIM:

To write a code in MATLAB to parse the NASA Thermodynamic data file and to calculate thermodynamic properties of various gas species.

OBJECTIVE:

• To write a function that extracts the 14 co-efficient of each species in data file.
• To calculate the enthalpy, entropy and specific heat for all species in data file.
• To calculate the molecular weight of each species and display it in command window.
• To plot the specific heat, Enthalpy and entropy for local temperature range specific for each species and save them.

MAIN PROGRAM:

% Code to parse the nasa thermodynamic data file and
%  to calculate thermodynamic properties of gas species.
clear all
close all
clc

% universal gas constant
R=8.134;

% Reading header
f1=fopen(\'THERMO.dat\',\'r\');
Header=fgetl(f1);

% Reading global temperature
Global_temp=fgetl(f1);
c=strsplit(Global_temp,\' \');
global_lowtemp=str2double(c{2});
global_midtemp=str2double(c{3});
global_maxtemp=str2double(c{4});

% Skipping 3 lines
fgetl(f1);
fgetl(f1);
fgetl(f1);

% creating for loop for reading each line of all species.

for i=1:53
    line=fgetl(f1);
    A=strsplit(line,\' \');
    species=A{1};
    Z=strfind(line,\'G\');
    local_lowtemp=line(Z+3:55);
    local_hightemp=line(Z+13:65);
    local_midtemp=line(Z+23:75);
    
    lt=str2double(local_lowtemp);
    ht=str2double(local_hightemp);
    mt=str2double(local_midtemp);
    
    Line_1=fgetl(f1);
    L1=strfind(Line_1,\'E\');
    c1=str2double(Line_1(1:(L1(1)+3)));
    c2=str2double(Line_1((L1(1)+4:(L1(2)+3))));
    c3=str2double(Line_1((L1(2)+4:(L1(3)+3))));
    c4=str2double(Line_1((L1(3)+4:(L1(4)+3))));
    c5=str2double(Line_1((L1(4)+4:(L1(5)+3))));
    
    Line_2=fgetl(f1);
    L2=strfind(Line_2,\'E\');
    c6=str2double(Line_2(1:(L2(1)+3)));
    c7=str2double(Line_2((L2(1)+4:(L1(2)+3))));
    c8=str2double(Line_2((L2(2)+4:(L1(3)+3))));
    c9=str2double(Line_2((L2(3)+4:(L1(4)+3))));
    c10=str2double(Line_2((L2(4)+4:(L1(5)+3))));
    
    Line_3=fgetl(f1);
    L3=strfind(Line_3,\'E\');
    c11=str2double(Line_3(1:(L3(1)+3)));
    c12=str2double(Line_3((L3(1)+4:(L1(2)+3))));
    c13=str2double(Line_3((L3(2)+4:(L1(3)+3))));
    c14=str2double(Line_3((L3(3)+4:(L1(4)+3))));
    
    % Ranging the temperature
    T=linspace(lt,mt,500);
    % specific heat
    cp=specific_heat(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,T,R,lt,ht,mt);
    % enthalpy
    H=enthalpy(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,T,R,lt,ht,mt);
    % entropy
    S=entropy(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,T,R,lt,ht,mt);
    % Creating folder to save plots of each species
     spec=char(species);
   fd=sprintf(\'%s\',spec);
   cur_dir=pwd;
   mkdir(spec);
   cd(spec);
   
   % plotting
   figure(1)
   plot( T,cp,\'linewidth\',2,\'color\',\'r\')
   xlabel(\'Temperature in K\')
   ylabel(\'specific heat in kJ\')
   title(species)
   
    figure(2)
   plot( T,S,\'linewidth\',2,\'color\',\'r\')
   xlabel(\'Temperature in K\')
   ylabel(\'Entropy in kJ/(mol k)\')
   title(species)
   
    figure(3)
   plot( T,H,\'linewidth\',2,\'color\',\'r\')
   xlabel(\'Temperature in K\')
   ylabel(\'Enthalpy in kJ/(mol)\')
   title(species)
   
   % saving plot in folder
   saveas(figure(1),\'cp.png\');
    saveas(figure(2),\'s.png\');
     saveas(figure(3),\'h.png\');
     cd(cur_dir);
     
   %Calculatimg Molecular weight
   h=0;c=0;n=0;o=0;ar=0;
   for k=1:length(species)
       if strcmp(species(k),\'H\')==1
           if k==length(species)
               h=h+1;
           else
               if isstrprop(species(k+1),\'digit\')==1
                   h=h+str2double(species(k+1));
               elseif isstrprop(species(k+1),\'digit\')==0
                   h=h+1;
               end
           end
       elseif strcmp(species(k),\'C\')==1
           if k==length(species)
           c=c+1;
       else
           if isstrprop(species(k+1),\'digit\')==1
               c=c+str2double(species(k+1));
           elseif isstrprop(species(k+1),\'digit\')==0
               c=c+1;
           end
           end
       elseif strcmp(species(k),\'N\')==1
           if k==length(species)
               n=n+1;
           
           else
               if isstrprop(species(k+1),\'digit\')==1
                   n=n+str2double(species(k+1));
               elseif isstrprop(species(k+1),\'digit\')==0
                   n=n+1;
               end
           end
   
           elseif strcmp(species(k),\'O\')==1
               if k==length(species)
                   o=o+1;
               else
                   if isstrprop(species(k+1),\'digit\')==1
                       o=o+str2double(species(k+1));
                   elseif isstrprop(species(k+1),\'digit\')==0
                       o=o+1;
                   end
               end
               elseif strcmp(species(k),\'AR\')==1
                   if k==length(species)
                       ar=ar+1;
                   else
                       if isstrprop(species(k+1),\'digit\')==1
                           ar=ar+str2double(species(k+1));
                       elseif isstrprop(species(k+1),\'digit\')==0
                           ar=ar+1;
                       end
                   end
                   elseif isstrprop(species(k),\'digit\')==1
                       continue;
  
       end
   end

H=1;C=12;N=14;O=16;AR=40;
molecular_weight=H*h+C*c+N*n+O*o+AR*ar;

fprintf(\'Molecularweightof %s:%f \\n\',species,molecular_weight);


end

FUNCTION PROGRAM FOR SPECIFIC HEAT (cp):

% Function to find specific heat
function cp=specific_heat(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,T,R,lt,ht,mt)
if T>=mt
    cp=(c1+c2.*(T)+c3.*(T.^2)+c4.*(T.^3)+c5.*(T.^4))*R;
else 
    cp=(c8+c9.*(T)+c10.*(T.^2)+c11.*(T.^3)+c12.*(T.^4))*R;
end

FUNCTION PROGRAM FOR ENTHALPY (H):

% Function to find enthalpy
function H=enthalpy(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,T,R,lt,ht,mt)
if T>=mt
    H=(c1.*T+c2.*(T.^2)./2+c3.*(T.^3)./3+c4.*(T.^4)./4+c5.*(T.^5)./5+c6)*R;
else
    H=(c8.*T+c9.*(T.^2)./2+c10.*(T.^3)./3+c11.*(T.^4)./4+c12.*(T.^5)./5+c13)*R;
end

FUNCTION PROGRAM FOR ENTROPY (S):

% Function to find entropy
function S=entropy(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,T,R,lt,ht,mt)
if T>=mt
    S=(c1.*log(T)+c2.*T+c3.*(T.^2)./2+c4.*(T.^3)./3+c5.*(T.^4)./4+c7)*R;
else
    S=(c8.*log(T)+c9.*T+c10.*(T.^2)./2+c11.*(T.^3)./3+c12.*(T.^4)./4+c14)*R;
end

EXPLANATION:

• The program uses the command fopen to open the thermodynamic file[thermo.dat] and stores it in file id called ‘f1’.
• ‘fgetl’ of f1 command takes the values of first line of file.

• Since temperatures are string weuse ‘strsplit’ command that splits the string into specified delimeter.
• Here ‘space’ is the delimiter,after splitting the values are assingned to variables.
• Using str2double command string value is converted into numeric value.
• The local temperatures are assigned to variables and linspace command used to take 500 values between local low to high temperatures.
• In the above code we are using ‘fgetl’command to extract coefficients.
• We are using strfind command to locate ‘E’ in each line.
• Strings are converted into numeric by str2double command
• Coefficient are assigned as c1,c2,c3,upto c14
• For loop is created to read the coefficient of all species.
• By creating separate function for specific heat(cp) , enthalpy(H) and entropy(S).
• We can calculate values of cp,H,S values by calling function.
• Command ‘cd(pwd)’ helps to open a present working directory.
• Mkdir command makes a new folder with name specified inside brackets.
• The molecular weight is also calculated by using if , for and elseif loop command.
• Finally the molecular weight of all species is displayed in command window.

FORMULA USED:

OUTPUT PLOTS FOR CO2,O2,N2:

CO2:

Specific heat:

Enthalpy:

Entropy:

O2:

Specific heat;

Enthalpy:

Entropy:

N2:

Specific heat:

Enthalpy:

Entropy:

SCREENSHOT OF FOLDERS SAVED :

MOLECULAR WEIGHT OF ALL SPECIES DISPLAYED: