Parsing Nasa Thermodynamic data using MATLAB

Objective : To write a program to read data from the given THERMO.dat file and to use this data to calculate specific heat , enthalpy, entropy and molecular weight if the species using the given formulae.

Procedure :

Step 1 : Initially to open the file fopen command is used.

fgetl command is used to read the line.

strfind command is used to find the position of 'E'.

Step 2 : Global temperature range is stored in the form of numbers using str2num command

Step 3 : Using for loop, data is read for each 53 species. For a given species, data is divided into 4 lines from line0 to line3.

line 0 contains the name of the species , the name of the species is extracted using strtok command

line1 to line3 contains 14 coefficients for the species out of which 7 are high temperature coefficients and remaining 7 are low temp coeffs.

These coeffs are extracted and are converted from string to numbers and then stored in variables ( l1, l2,........,l7 and h1,h2,......h7)

Step 4 : Afte extracting the coeffs they are used in the given formule to calculate properties and molecular weight.

while defining function to calculate molecular weight, strcmpi command is used instead of strcmp since strcmp command compares the strings irrespective of the case of string.

Step 5 : New directory is used by making use of mkdir command.

Formula :

Definitions :

Enthalpy is a thermodynamic property of a system. It is the sum of the internal energy added to the product of the pressure and volume of the system.

Entropy the measure of a system's thermal energy per unit temperature that is unavailable for doing useful work. Because work is obtained from ordered molecular motion, the amount of entropy is also a measure of the molecular disorder, or randomness, of a system.

specific heat is the amount of heat per unit mass required to raise the temperature by one degree Celsius.

MATLAB CODE :

close all
clear all
clc


R = 8.314;
f1 = fopen('THERMO.txt','r')

file_header = fgetl(f1)

temp_range = fgetl(f1)

Global_temp = strsplit(temp_range,' ')
Global_low_temp = str2double(Global_temp{1})
Global_mid_temp = str2double(Global_temp{2})
Global_high_temp = str2double(Global_temp{3})

T= linspace(Global_low_temp,Global_high_temp,500)

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


for i = 1 : 53 
    
line0 = fgetl(f1)             % first line of data 

species_name = strtok(line0)     


% reading the first three lines which contains 7 high temperature coeffs and 7 low temperature coeffs
line1 = fgetl(f1)
line2 = fgetl(f1)
line3 = fgetl(f1)


% Finding the position of E
A = strfind(line1, 'E' )
B = strfind(line2, 'E' )
C = strfind(line3, 'E' )


    % Extracting Higher temperature coeffs
   
    l1 = str2double(line1(1:A(1)+3));
    l2 = str2double(line1(A(1)+4:A(2)+3));
    l3 = str2double(line1(A(2)+4:A(3)+3));
    l4 = str2double(line1(A(3)+4:A(4)+3));
    l5 = str2double(line1(A(4)+4:A(5)+3));
    l6 = str2double(line2(1:B(1)+3));
    l7 = str2double(line2(B(1)+4:B(2)+3));
    
    % extracting lower temperature coeffs
    
    h1 = str2double(line2(B(2)+4:B(3)+3));
    h2 = str2double(line2(B(3)+4:B(4)+3));
    h3 = str2double(line2(B(4)+4:B(5)+3));
    h4 = str2double(line3(1:C(1)+3));
    h5 = str2double(line3(C(1)+4:C(2)+3));
    h6 = str2double(line3(C(2)+4:C(3)+3));
    h7 = str2double(line3(C(3)+4:C(4)+3));
    
    
    
    cp = specific_heat(l1,l2,l3,l4,l5,h1,h2,h3,h4,h5,T,R, Global_mid_temp);
    s= entropy(l1,l2,l3,l4,l5,l7,h1,h2,h3,h4,h5,h7,T,R,Global_mid_temp);
    h= enthalpy(l1,l2,l3,l4,l5,l6,h1,h2,h3,h4,h5,h6,T,R,Global_mid_temp);
    molecular_weigth(i) = molecular_wt(species_name)
    plotting(cp,h,s,T,species_name)
    
    
end

code for specific heat :

function cp = specific_heat(l1,l2,l3,l4,l5,h1,h2,h3,h4,h5,T,R, Global_mid_temp)

if T<=Global_mid_temp
    
    cp = R*(l1 + l2.*(T) + l3.*(T.^2) + l4.*(T.^3) + l5.*(T.^4));
    
else
    
    cp = R*(h1 + h2.*(T) + h3.*(T.^2) + h4.*(T.^3) + h5.*(T.^4) );
    
end

code for enthalpy :

function h= enthalpy(l1,l2,l3,l4,l5,l6,h1,h2,h3,h4,h5,h6,T,R,Global_mid_temp)

if T<= Global_mid_temp
    
    h = R*(l1.*T + l2.*(T.^2)./2 + l3.*(T.^3)./3 + l4.*(t.^4)./4  + l5.*(t.^5)./5 + l6 );
    
else
    
    h = R*( h1.*T + h2.*(T.^2)./2  + h3.*(T.^3)./3  + h4.*(T.^4)./4 + h5.*(T.*5)./5 + h6 )
    
end

code for entropy :

function s= entropy(l1,l2,l3,l4,l5,l7,h1,h2,h3,h4,h5,h7,T,R,Global_mid_temp)

if T<=Global_mid_temp
    
    s = R*(l1.*log(T)+ l2.*T + l3.*(T.^2)./2 + l4.*(T.^3)./3 + l5.*(T.^4)./4 + l7 );
    
else
    
    s= R*( h1.*log(T)+ h2.*T + h3.*(T.^2)./2 + h4.*(T.^3) ./3 + h5.*(T.^4)./4 + h7 );
    
end

code for molecular weight :

function [MW]=molecular_wt(species_name)

Element = { 'H' 'C' 'N' 'O' 'R' } ;

fprintf(' molecular weight of ')
fprintf('species_name')
fprintf(' : ')

atomic_wt = [ 1 12 14 16 40 ]    % atomic weight of each element

MW = 0

for i = 1:length(species_name)
for j = 1:length(Element)
    
    if strcmpi(species_name(i),Element(j))
        MW = MW + atomic_wt(j)
        position=j                   % elements stored in position
    end
end
    
    n = str2double(species_name(i))
    
    if n>1
        
        MW = MW + atomic_wt(position)*(n-1);
    end
end
        fprintf('%d',MW)
        disp(' ')
       
end

plots :

1. Enthalpy vs Temperature :

2. Entrophy vs Temperature

3. Specific heat vs Temperature

Reference :

https://in.mathworks.com/help/matlab/ref/strcmp.html

https://in.mathworks.com/matlabcentral/answers/86318-separating-values-in-cell-arrays