Challenge

Project 1 - Parsing NASA thermodynamic data

SSikkandar Faiz18 Mar 2022

THINGS TO FIND:

1.read the thermo file of nasa to find the contents of it.

2.once after figuring out the data and how to do the code take out the important contents of the file and try to come up with a logic to bring out function

3.take out the coefficients of the given gas species to find its specific heat(cp), entropy(S), enthalpy(H) for all the gas elements.

plot the graph for all of them for their cp,H,S vs temperature.

4.files containing the plots of all gas elements are to be saved in a new folder.

POLYNOMIAL FUNCTOION FOR FINDING THE cp,H AND S are:

................................................................................................................................................................................................................................................................................................................................................................

MAIN CODE:

clear all 
close all 
clc
X = fopen('THERMO.dat','r');
%HEADING
line1 = fgetl(X);
%global temperature value
line2 = fgetl(X);

%NEXT THREE LINES COMMENTS
line3 = fgetl(X);
line4 = fgetl(X);
line5 = fgetl(X);

%getting global temperature value into numerical format
gtv = strsplit(line2,' ');
global_low_temp = str2num(gtv{2});
global_mid_temp = str2num(gtv{3});
global_hig_temp = str2num(gtv{4});
%temperature(T) and R value(universal gas constant)
T = linspace(global_low_temp,global_hig_temp,500);
R = 8.314 


%getting into species and their coefficients and making the code to run 53
%times as we have 53 species in the data
for i = 1:53
    %to get next line we use fgetl which here is the o, oxygen.
    species = fgetl(X);
    %we can use strtok command to take out the name of the species
    element_name = strtok(species)
    %next lines contains the coefficients of the species needed to know the
    %cp, H AND S
    coeff1 = fgetl(X);
    A = strfind(coeff1,'E')
    coeff2 = fgetl(X);
    B = strfind(coeff2,'E')
    coeff3 = fgetl(X);
    C = strfind(coeff3,'E')
    %now we are to use a and b to identify the first seven and the second
    %seven coeffs
    a1 = str2num(coeff1(1:A(1)+3));
    a2 = str2num(coeff1(A(1)+4:A(2)+3));
    a3 = str2num(coeff1(A(2)+4:A(3)+3));
    a4 = str2num(coeff1(A(3)+4:A(4)+3));
    a5 = str2num(coeff1(A(4)+4:A(5)+3)); 
    
    a6 = str2num(coeff2(1:B(1)+3));
    a7 = str2num(coeff2(B(1)+4:B(2)+3));

    b1 = str2num(coeff2(B(2)+4:B(3)+3));
    b2 = str2num(coeff2(B(3)+4:B(4)+3));
    b3 = str2num(coeff2(B(4)+4:B(5)+3));

    b4 = str2num(coeff3(1:C(1)+3));
    b5 = str2num(coeff3(C(1)+4:C(2)+3));
    b6 = str2num(coeff3(C(1)+4:C(2)+3));
    b7 = str2num(coeff3(C(1)+4:C(2)+3));

    %calcuating specific heat
    cp = specific_heat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,global_mid_temp);
    %calculating entropy
    S = entropy(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T,R,global_mid_temp);
    %calculation enthalpy
    H = enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T,R,global_mid_temp);

    %molecular weight calculation for the elements
    MW = molecular_weight(element_name);

    %save plots in different folders
    cd('C:\Users\SIKKANDAR FAIZ\Documents\MATLAB\THERMODDYNAMIC DATA')
    dir_id = element_name;
    mkdir(dir_id)

    figure(1)
    plot(T,cp,'LineWidth',3,'color','r')
    title(sprintf('specific heat vs temperature of %s',element_name))
    xlabel('temperature')
    ylabel('specific heat')
    grid on


    figure(2)
    plot(T,S,'LineWidth',3,'color','b')
    title(sprintf('entropy vs temperature of %s',element_name))
    xlabel('temperature')
    ylabel('entropy')
    grid on
    

    figure(3)
    plot(T,H,'LineWidth',3,'color','g')
    title(sprintf('enthalpy vs temperature of %s',element_name))
    xlabel('temperature')
    ylabel('enthalpy')
    grid on
    
    %pause

    pause(0.0001)
    
    cd(dir_id);
    saveas(1,'specific heat.png');
    saveas(2,'entropy.png');
    saveas(3,'enthalpy.png');
end
fclose(X);

MAIN CODE EXPLANATION:

i will try to explain every single line as i can:

first like every time clear all, close all, clc to refresh the system for the upcoming code to function effectively.

clear all

close all

clc

once after the above, we move on to taking out the data from the given thermodata. which contains the coefficients of the gas species to find the specific heat, entropy and enthalpy of all them by using the nasa polynomials. X = fopen('THERMO.dat','r') this opens the given file data for further processing.

X = fopen('THERMO.dat','r');

next we use the fgetl command to take out data from next line. we are assigning the second line we get to a name called line1.this line contains the name of file in the given data.

%HEADING

line1 = fgetl(X);

we are to use the same command to take out the next line, which here contains the global temperatue value. it is assigned name line2

%global temperature value

line2 = fgetl(X);

same command for the next three lines which contains only comments, which aren't much use to the code. assigned some line names.

%NEXT THREE LINES COMMENTS

line3 = fgetl(X);

line4 = fgetl(X);

line5 = fgetl(X);

we now need to unpack the values of global temperature as they are in an array format and are reoresented as strings. we are using the strsplit command, which split the data of global temperature where it finds the space(' ') as we made it to perform.

%getting global temperature value into numerical format

gtv = strsplit(line2,' ');

after that we will still get the values as strings but in a split way. so we are to make all the strings we got into numbers. we are using str2num command which will help us in what we want. after which we are assigning the names of the three got number values(300,1000 and 5000) as global_low_temp, global_mid_temp and global_hig_temp

global_low_temp = str2num(gtv{2});

global_mid_temp = str2num(gtv{3});

global_hig_temp = str2num(gtv{4});

the temperature we got is then arranged as we need it, by using linspace command. and R value which is the universal gas constant is took for consideration.

%temperature(T) and R value(universal gas constant)

T = linspace(global_low_temp,global_hig_temp,500);

R = 8.314

we are now to get to the main part of the code. we need to find cp,h and S of all the gas elements present in the given data. when we count the gas elements we will be seeing there are 53 elements. we are hence making a code to find the values of all the elements. we are to run the code 53 time to find all values. function start.

%getting into species and their coefficients and making the code to run 53

%times as we have 53 species in the data

for i = 1:53

we then come up with the same fgetl command to get to next line. which here contains the name of the element(gas)

%to get next line we use fgetl which here is the o, oxygen.

species = fgetl(X);

we then use the strtok command to take out the name of the element.

%we can use strtok command to take out the name of the species

element_name = strtok(species)

we then move onto the next three lines which contains the data for the calculations using fgetl. once getting all the values in the next line of the element name, we are to unpack it as they are packed together. here we use strfind command to id a letter that is same to all the coefficients.we are to do the same to the next two line which contain the coefficients. which here we find is'E'. we are to assign the index position of all the E in the three lines.the E in three lines are differentiated with the help of vae=riables A,B and C.

%next lines contains the coefficients of the species needed to know the

%cp, H AND S

coeff1 = fgetl(X);

A = strfind(coeff1,'E')

coeff2 = fgetl(X);

B = strfind(coeff2,'E')

coeff3 = fgetl(X);

C = strfind(coeff3,'E')

with the help of the index position of the E we then splice the 14 values that are the coefficients of the gas element to find it cp,H,S. first seven represents the high temp coefficients and the next seven are low. we differentiate them as a and b.

%now we are to use a and b to identify the first seven and the second

%seven coeffs

a1 = str2num(coeff1(1:A(1)+3));

a2 = str2num(coeff1(A(1)+4:A(2)+3));

a3 = str2num(coeff1(A(2)+4:A(3)+3));

a4 = str2num(coeff1(A(3)+4:A(4)+3));

a5 = str2num(coeff1(A(4)+4:A(5)+3));

a6 = str2num(coeff2(1:B(1)+3));

a7 = str2num(coeff2(B(1)+4:B(2)+3));

b1 = str2num(coeff2(B(2)+4:B(3)+3));

b2 = str2num(coeff2(B(3)+4:B(4)+3));

b3 = str2num(coeff2(B(4)+4:B(5)+3));

b4 = str2num(coeff3(1:C(1)+3));

b5 = str2num(coeff3(C(1)+4:C(2)+3));

b6 = str2num(coeff3(C(1)+4:C(2)+3));

b7 = str2num(coeff3(C(1)+4:C(2)+3));

we are using the specific heat sub code to calculate the cp value of gases. will explain it below.

%calcuating specific heat

cp = specific_heat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,global_mid_temp);

we are using the entropy sub code to calculate the S value of gases. will explain it below.

%calculating entropy

S = entropy(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T,R,global_mid_temp);

we are using the enthalpy sub code to calculate the H value of gases. will explain it below.

%calculation enthalpy

H = enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T,R,global_mid_temp);

we are using the molecular weight sub code to calculate the MW value of gases. will explain it below.

%molecular weight calculation for the elements

MW = molecular_weight(element_name);

the plots required are saved in a new folder.

%save plots in different folders

cd('C:\Users\SIKKANDAR FAIZ\Documents\MATLAB\THERMODDYNAMIC DATA')

dir_id = element_name;

mkdir(dir_id)

plot one contains the cp vs T for the gas elements

figure(1)

plot(T,cp,'LineWidth',3,'color','r')

title(sprintf('specific heat vs temperature of %s',element_name))

xlabel('temperature')

ylabel('specific heat')

grid on

plot two contains the S vs T for the gas elements

figure(2)

plot(T,S,'LineWidth',3,'color','b')

title(sprintf('entropy vs temperature of %s',element_name))

xlabel('temperature')

ylabel('entropy')

grid on

plot three contains the H vs T for the gas elements

figure(3)

plot(T,H,'LineWidth',3,'color','g')

title(sprintf('enthalpy vs temperature of %s',element_name))

xlabel('temperature')

ylabel('enthalpy')

grid on

we are pausing to have a better view of the plots

%pause

pause(0.0001)

the names if the plots are given in order to be saved in their own folder.

cd(dir_id);

saveas(1,'specific heat.png');

saveas(2,'entropy.png');

saveas(3,'enthalpy.png');

function end

end

close file

fclose(X);

CODE FOR MOLECULAR WEIGHT CALCULATION:

function MW = molecular_weight(element_name)
main_atoms = ['O','C','N','H','A'];
fprintf('MOLECULAR WEIGHT OF ')
fprintf(element_name);
fprintf(' : ');
atomic_weight = [16 12 14 1 40];
MW = 0;
for i = 1:length(element_name);
    for j = 1:length(main_atoms);
        if strcmp(element_name(i),main_atoms(j))
            MW = MW+atomic_weight(j);
            position = j;
        end
    end
    n = str2double(element_name(i));
    if n>1
        MW = MW + atomic_weight(position)*(n-1);
    end
end
fprintf('%d',MW);
disp('');
end

MOLECULAR WEIGHT CALCULATION CODE EXPLANATION:

we start with the function command

function MW = molecular_weight(element_name)

we consider all the required main atoms

main_atoms = ['O','C','N','H','A'];

we then name how our output should be displayed

fprintf('MOLECULAR WEIGHT OF ')

fprintf(element_name);

fprintf(' : ');

the atomic weight of the main elements are taken into consideration as the molecular weight formula requires them

atomic_weight = [16 12 14 1 40];

MW first is given value 0

MW = 0;

for i go upto 53, for every 53 moves of i go verify all 5 main elements of j. in which if the code finds the name of element match main atom then add mw and the atomic weight of matched main atom and position it back to j.

for i = 1:length(element_name);

for j = 1:length(main_atoms);

if strcmp(element_name(i),main_atoms(j))

MW = MW+atomic_weight(j);

position = j;

end

n is the number molecules of the element

n = str2double(element_name(i));

if n>1

MW = MW + atomic_weight(position)*(n-1);

end

fprintf('%d',MW);

disp('');

end

CODE FOR SPECIFIC HEAT:

function cp = specific_heat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,global_mid_temp)
%comparing with global_mid_temp 
if T > global_mid_temp
    cp = R.*(a1 + a2.*T + a3.*T.^2 + a4.*T.^3 + a5.*T.^4 );
else
    cp = R.*(b1 + b2.*T + b3.*T.^2 + b4.*T.^3 + b5.*T.^4 );
end

EXPLANATION:

from the above given polynomials for specific heat, the taken out coefficients of the gas element along with the temperature, universal gas constant and global mid temp, the code is made to find the cp for the gases.if the temperature(T) is below global mid temp then the values of b which are the low temperature coefficients are taken into consideration for code, if not then a values which are high temp coefficients are considered.

CODE FOR ENTROPY:

function S = entropy(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T,R,global_mid_temp)
%comparing T with global_mid_temp
if T > global_mid_temp
    S = R.*(a1.*log(T) + a2.*T + a3.*(T.^2./2) + a4.*(T.^3./3) + a5.*(T.^4./4) + a7);
else
    S = R.*(b1.*log(T) + b2.*T + b3.*(T.^2./2) + b4.*(T.^3./3) + b5.*(T.^4./4) + b7);
end

EXPLANATION:

from the above given polynomials for entropy, the taken out coefficients of the gas element along with the temperature, universal gas constant and global mid temp, the code is made to find the S for the gases.if the temperature(T) is below global mid temp then the values of b which are the low temperature coefficients are taken into consideration for code, if not then a values which are high temp coefficients are considered.

CODE FOR ENTHALPY:

function H = enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T,R,global_mid_temp)
%comparing with global_mid_temp
if T > global_mid_temp
    H = R.*T.*(a1 + a2.*(T./2) + a3.*(T.^2./3) + a4.*(T.^3./4) + a5.*(T.^4./5) + a6./T);
else
    H = R.*T.*(b1 + b2.*(T./2) + b3.*(T.^2./3) + b4.*(T.^3./4) + b5.*(T.^4./5) + b6./T);
end

EXPLANATION:

from the above given polynomials for enthalpy, the taken out coefficients of the gas element along with the temperature, universal gas constant and global mid temp, the code is made to find the H for the gases.if the temperature(T) is below global mid temp then the values of b which are the low temperature coefficients are taken into consideration for code, if not then a values which are high temp coefficients are considered.

PLOTS OF O2,N2 AND CO2:

PLOTS OF O2:

................................................................................................................................................................................

PLOTS OF N2:

................................................................................................................................................................................

PLOTS OF CO2:

MOLECULAR WEIGHT OF ELEMENTS:

MOLECULAR WEIGHT OF O : 16
MOLECULAR WEIGHT OF O2 : 32
MOLECULAR WEIGHT OF H : 1
MOLECULAR WEIGHT OF H2 : 2
MOLECULAR WEIGHT OF OH : 17
MOLECULAR WEIGHT OF H2O : 18
MOLECULAR WEIGHT OF HO2 : 33
MOLECULAR WEIGHT OF H2O2 : 34
MOLECULAR WEIGHT OF C : 12
MOLECULAR WEIGHT OF CH : 13
MOLECULAR WEIGHT OF CH2 : 14
MOLECULAR WEIGHT OF CH2(S) : 14
MOLECULAR WEIGHT OF CH3 : 15
MOLECULAR WEIGHT OF CH4 : 16
MOLECULAR WEIGHT OF CO : 28
MOLECULAR WEIGHT OF CO2 : 44
MOLECULAR WEIGHT OF HCO : 29
MOLECULAR WEIGHT OF CH2O : 30
MOLECULAR WEIGHT OF CH2OH : 31
MOLECULAR WEIGHT OF CH3O : 31
MOLECULAR WEIGHT OF CH3OH : 32
MOLECULAR WEIGHT OF C2H : 25
MOLECULAR WEIGHT OF C2H2 : 26
MOLECULAR WEIGHT OF C2H3 : 27
MOLECULAR WEIGHT OF C2H4 : 28
MOLECULAR WEIGHT OF C2H5 : 29
MOLECULAR WEIGHT OF C2H6 : 30
MOLECULAR WEIGHT OF CH2CO : 42
MOLECULAR WEIGHT OF HCCO : 41
MOLECULAR WEIGHT OF HCCOH : 42
MOLECULAR WEIGHT OF H2CN : 28
MOLECULAR WEIGHT OF HCN : 27
MOLECULAR WEIGHT OF HNO : 31
MOLECULAR WEIGHT OF N : 14
MOLECULAR WEIGHT OF NNH : 29
MOLECULAR WEIGHT OF N2O : 44
MOLECULAR WEIGHT OF NH : 15
MOLECULAR WEIGHT OF NH2 : 16
MOLECULAR WEIGHT OF NH3 : 17
MOLECULAR WEIGHT OF NO : 30
MOLECULAR WEIGHT OF NO2 : 46
MOLECULAR WEIGHT OF HCNO : 43
MOLECULAR WEIGHT OF HOCN : 43
MOLECULAR WEIGHT OF HNCO : 43
MOLECULAR WEIGHT OF NCO : 42
MOLECULAR WEIGHT OF CN : 26
MOLECULAR WEIGHT OF HCNN : 41
MOLECULAR WEIGHT OF N2 : 28
MOLECULAR WEIGHT OF AR : 40
MOLECULAR WEIGHT OF C3H8 : 44
MOLECULAR WEIGHT OF C3H7 : 43
MOLECULAR WEIGHT OF CH3CHO : 44
MOLECULAR WEIGHT OF CH2CHO : 43

FOLDER OF PLOTS:

ERRORS:

there was some minor file name issues. sorted it later.

CONCLUSION:

With help from previous class videos and many youtube videos and official matlab website, the code ran well and the cp,H, S are found with the help of the thermodata.

Project 1 - Parsing NASA thermodynamic data

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.