Project 1 - Parsing NASA thermodynamic data

Aim - Write a code in Matlab to parse the NASA thermodynamic data file and then calculate thermodynamic properties of various gas species.

The file from which the data is taken is THERMO.dat

The formulas used for finding the Specific Heat is -

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

Specific Heat Calculation

%specific heat (cp) calculation function

function cp = specificheat_calc(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,R,tempmid)

if temp <=  tempmid
    
    cp = R.*(a8 + (a9.*temp) + (a10.*(temp).^2) + (a11.*(temp).^3) + (a12.*(temp).^4));
else
    cp = R.*(a1 + (a2.*temp) + (a3.*(temp).^2) + (a4.*(temp).^3) + (a5.*(temp).^4));
end

The formulas used for finding the Entropy is -

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

Entropy Calculation

%Entropy (S) calculation function

function S = Entropy_calc(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,R,tempmid)

if temp <=  tempmid
    
    S = R.*(a8.*log(temp) + (a9.*temp) + (a10.*((temp).^2)/2) + (a11.*((temp).^3)/3) + (a12.*((temp).^4)/4) + a14);
else
    S = R.*(a1.*log(temp) + (a2.*temp) + (a3.*((temp).^2)/2) + (a4.*((temp).^3)/3) + (a5.*((temp).^4)/4) + a7);
end

The formulas used for finding the Enthalpy is -

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

Enthalpy Calculation

%Enthalpy (H) calculation function

function H = Enthalpy_calc(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,temp,R,tempmid)

if temp <=  tempmid
    
    H = R.*temp.*(a8 + (a9.*(temp)/2) + ((a10.*(temp).^2)/3) + (a11.*((temp).^3)/4) + (a12.*((temp).^4)/5) + ((a13./temp)));
else
    H = R.*temp.*(a1 + (a2.*(temp)/2) + ((a3.*(temp).^2)/3) + (a4.*((temp).^3)/4) + (a5.*((temp).^4)/5) + ((a6./temp)));
end

Molecular Weight Calculation

%Molecular Weight Calculation Function

function W = Molecular_weight(Species_SP)

atomic_name  = ['H','C','O','N','Ar'];

% fprintf('Molecular Weight of ');
% fprintf(Species_SP);
% fprintf(' : ');

Atomic_Weight = [1.008 12.011 15.999 14.007 39,948];

W = 0;

for i = 1:length(Species_SP)
    for j = 1:length(atomic_name)
        if strcmp(Species_SP(i),atomic_name(j))
            W = W + Atomic_Weight(j);
            position = j;
        end
    end
    n = str2num(Species_SP);
    
    if n>1
        W = W + Atomic_Weight(position)*(n-1);
    end
end
    fprintf('molecular weight %s : ',Species_SP);
    fprintf('%f',W);
    disp(' ');
end

Main Code

clear all
close all
clc

%universal gas constant J/(mol-K)
    R = 8.314;
    
%making directory
mkdir 'D:SoftwaresMATLAB'

%reading data file - Nasa Thermo Data
f1 = fopen('THERMO.dat','r');
gt = fgetl(f1);

%obtaining temps from reading header information
temp = fgetl(f1);

%splitting
t = strsplit(temp);

%converting strings to numbers
templow = str2num(t{2});
tempmid = str2num(t{3});
temphigh = str2num(t{4});

%temp renge and gas constant
temp = linspace(templow,temphigh,1000);

    
%skipping comment lines in file
for i = 1:3
    gt = fgetl(f1);
end

%reading Species
for i = 1:53
    gt = fgetl(f1);
    A = strsplit(gt,' ');
    Species_SP = (A{1});
    
    line_1 = fgetl(f1);
    a = strfind(line_1,'E');
    
    %converting string to number
    lowtemp = str2num(A{end-3});
    midtemp = str2num(A{end-2});
    hightemp = str2num(A{end-1});
    
    %extracting high and low temp coefficients
    
    %Extracting first 7 (high) temp coeffs
    
    a1 = line_1(1:a(1) + 3);
    a1 = str2num(a1);
    
    a2 = line_1(a(1) + 4:a(2) + 3);
    a2 = str2num(a2);
    
    a3 = line_1(a(2) + 4:a(3) + 3);
    a3 = str2num(a3);
    
    a4 = line_1(a(3) + 4:a(4) + 3);
    a4 = str2num(a4);
    
    a5 = line_1(a(4) + 4:a(5) + 3);
    a5 = str2num(a5);
    
    line_2 = fgetl(f1);
    b = strfind(line_1,'E');
    
    a6 = line_2(1:b(1) + 3);
    a6 = str2num(a6);
    
    a7 = line_2(a(1) + 4:b(2) + 3);
    a7 = str2num(a7);
    
    %extracting second 7 (low) temp coeffs
    
    a8 = line_2(b(2) + 4:b(3) + 3);
    a8 = str2num(a8);
    
    a9 = line_2(b(3) + 4:b(4) + 3);
    a9 = str2num(a9);
    
    a10 = line_2(b(4) + 4:b(5) + 3);
    a10 = str2num(a10);
    
    line_3 = fgetl(f1);
    c = strfind(line_1,'E');
    
    a11 =  line_3(1:c(1) + 3);
    a11 = str2num(a11);
   
    a12 = line_3(c(1) + 4:c(2) + 3);
    a12 = str2num(a12);
    
    a13 = line_3(c(2) + 4:c(3) + 3);
    a13 = str2num(a13);
    
    a14 = line_3(c(3) + 4:c(4) + 3);
    a14 = str2num(a14);
    
    %specific heat calculation
    cp = specificheat_calc(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,R,tempmid);
    
    %entropy calculation
    S = Entropy_calc(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,R,tempmid);
    
    %enthalpy calculation
    H = Enthalpy_calc(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,temp,R,tempmid);
    
    %calculating molecular weight of species
    W = Molecular_weight(Species_SP);
    
    %read and write molecular list
    Molecular_list = fopen('Molecular weight .txt','w');
    fprintf(Molecular_list,'Molecular weight of %s is %d ',Species_SP,W);
    fclose(Molecular_list);
    
    %plotting(CH,H,S,temp,Species_SP)
    %saving as a folder
    
    cur_dir = pwd;
    mkdir(Species_SP);
    cd(Species_SP);
    
    % plot 1 - Specific Heat (cp) vs Temperature Range (temp)
    
    figure(1)
    plot(temp,cp,'linewidth',2,'color','r');
    xlabel('Temp [K]');
    ylabel('Specific Heat [kJ]');
    title(sprintf('varying specific heat with temperature range of %s',Species_SP));
    grid on
    saveas(figure(1),'Specific Heat.jpg');
    
    %plot 2 - Enthalpy (H) vs Temperature Range (temp)
    
    figure(2)
    plot(temp,H,'linewidth',2,'color','g');
    xlabel('Temp [K]');
    ylabel('Enthalpy [kJ/mol]');
    title(sprintf('varying enthalpy with temperature range of %s',Species_SP));
    grid on
    saveas(figure(2),'Enthalpy.jpg');
    
    %plot 3 - Entropy (S) vs Temperature Range (temp)
    
    figure(3)
    plot(temp,S,'linewidth',2,'color','b');
    xlabel('Temp [K]');
    ylabel('Entropy [kJ/mol-K]');
    title(sprintf('varying entropy with temperature range of %s',Species_SP));
    grid on
    saveas(figure(3),'Entropy.jpg');
    
    cd(cur_dir);
    
    %Molecular Weight
    
    W = Molecular_weight(Species_SP)
    f2 = fopen('molecular mass','w');
    fprintf('%s n %s n %f n','species','Molecular Mass',W)
    fclose(f2);
end

Steps -

• First we input the parameter R and open the THERMO.dat file and read the file using 'fopen' command.
• We read the first line of the file using 'fgetl' command, read the temp file and split the temps using 'strsplit' command to convert it into a cell array.
• Then we convert this string into numbers using the 'str2num' command to find the low high and mid temps.
• Then assign temp range using the 'linspace' command. we run the loop to skip the command ines in the date file using 'fgetl' command in the loop.
• We run the loop to find the 53 elements then splitting into the cell array to find the species line then the species name, then the number of temps in the species line.
• When we read the new line.
• In order to read the coefficient in the data file we wil fidn the letter 'E' by using 'strfind' command to find before and after numbers of the coefficients.
• Like this we find the remaining 14 coefficients and convert them to numbers from strings.
• The first 7 coefficients are for Higher Temps, while the second 7 are for Lower Temps.
• Now we specify Cp, S, H, from the coefficients, then the function Molecular_weight will find the weights of the different species.
• We write the Molecular weight.txt ffile and print all the weights of the species using 'fprintf' and then close the txt file using 'fclose' command.
• We change the directory using 'cd' command and creating directory folder where we have to save the files using 'mkdir' command.
• Plot the Cp, S and H vs Temps for all species.
• Using 'saveas' ccommand to save the plots of different species in the directory folder.

Plots of O2-

Plots of N2-

Plots of CO2-

The command window values has been pasted word to word in a notepad which is given in the link - https://drive.google.com/file/d/1p05PbYDXpTgP6s5A61dyPvk-YU0q8lTl/view?usp=sharing

https://drive.google.com/drive/folders/19p6yQ1nhNgBlMzceebqv8eHTTeC6kA16?usp=sharing