Project 1 - Parsing NASA thermodynamic data

OBJECTIVE : 

Write a function that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the data file. 

 Here R is the universal gas constant, T is the temeprature

The FIRST 7 coefficients are HIGH-temperature coefficients and the SECOND 7 coefficients are LOW-temperature coefficients.

Also use the LOCAL TEMPERATURE RANGES which is unique for each individual species, not the GLOBAL TEMPERATURE RANGE.

1.  Calculate the molecular weight of each species and display it in the command window.
2. Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species.
3. Save the plots as images with appropriate names and dump them in separate folders for each species with suitable name. All these should be generated automatically.

THEORY & MATLAB CODE :

 Let us primarily understand what is enthalpy , entropy and specific heat.

Thermodynamics is basically the study of effects of heat, work and energy on the whole. Under thermodynamics there are different parameters that branch ,

Enthalpy - It is the main factor of thermodynamics. It is the thermodynamic quantity equivalent to the total heat content of a system. And this means if the energy is added, the enthalpy increases. If the energy is given off, then the enthalpy of the system decreases.

We represent enthalpy using the formulae ,

H = U + PV

• U = It is basically the heat content or the heat excreted by atoms and molecules at a more discrete level or a quantum level.
• P = Pressure
• V = Volume

Enthalpy basically is said to be the total amount of energy a certain object is said to possess based on their molecular attraction an energy is exhibited which is potential energy. And this is expressed as Internal Energy and pressure + Volume , energy used to establish a system from empty space.

Entropy - The amount of entropy is also a measure of the molecular disorder, or randomness, of a system. 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.

Enthalpy is the sum total of all the energies, whereas entropy is the measure of the change in enthalpy/temperature. 

Specific Heat - The specific heat is the amount of heat per unit mass required to raise the temperature by one degree Celsius. The relationship between heat and temperature change is usually expressed in the form shown below where c is the specific heat. The relationship does not apply if a phase change is encountered, because the heat added or removed during a phase change does not change the temperature.

Let us understand the code,

• Primarily we open the 'thermo.dat' file using the Fopen command and we assign it to a variable called 'nasa_data'.
• Now when we use fgetl , we get the first value of the 'thermo.dat' file., which is the title of the document.
• Now in order to obtain the Global temperature range, we get assign a variable called 'global_temp_range' and obtain the next field value in the document.
• But we are trying to split one string value into separate discrete string values using the command strsplit. When we use the strsplit, we obtain an array of values within the variable temp_scale. As a result, we can address the particular field value by indexing the value using the variable hande.
• We store the values of global_low_temp ,  global_med_temp & global_high_temp values inside each variable by indexing the variable handle and the values stored inside the array. And also we save the value inside the corresponding variables as integer values using 'str2num' function.
• Inorder to read out the comments , we create a for loop and read out the next three lines under the comments variable handle.

clear all
close all
clc

% Title
nasa_data = fopen('THERMO.dat');
Title = fgetl(nasa_data)

% Global Temperature Range
global_temp_range = fgetl(nasa_data);
temp_scale = strsplit(global_temp_range,' ')
global_low_temp = str2num(temp_scale{2})
global_med_temp = str2num(temp_scale{3})
global_high_temp = str2num(temp_scale{4})

% Comments 
for headers = 1 : 3
    comments = fgetl(nasa_data)
end


% O - Oxygen Coefficients for High Temperature & Low Temperature

% Local Temperature Ranges
for i = 1:53
    line_1 = fgetl(nasa_data)
    a = strsplit(line_1,' ')
    species=a{1}

    min_temp = str2double(a(length(a) - 3))
    max_temp = str2double(a(length(a) - 2))
    mid_temp = str2double(a(length(a) - 1))

    line_2 = fgetl(nasa_data)
    a1 = strfind(line_2,'E')
    h_t1 = str2num(line_2(1 : a1(1)+3))
    h_t2 = str2num(line_2(a1(1)+4 : a1(2)+3))
    h_t3 = str2num(line_2(a1(2)+4 : a1(3)+3))
    h_t4 = str2num(line_2(a1(3)+4 : a1(4)+3))
    h_t5 = str2num(line_2(a1(4)+4 : a1(5)+3))

    line_3 = fgetl(nasa_data)
    a2 = strfind(line_3,'E')
    h_t6 = str2num(line_3(1 : a2(1)+3))
    h_t7 = str2num(line_3(a2(1)+4 : a2(2)+3))
    l_t1 = str2num(line_3(a2(2)+4 : a2(3)+3))
    l_t2 = str2num(line_3(a2(3)+3 : a2(4)+3))
    l_t3 = str2num(line_3(a2(4)+4 : a2(5)+3))

    line_4 = fgetl(nasa_data)
    a3 = strfind(line_4,'E')
    l_t4 = str2num(line_4(1 : a3(1)+3))
    l_t5 = str2num(line_4(a3(1)+4 : a3(2)+3))
    l_t6 = str2num(line_4(a3(2)+4 : a3(3)+3))
    l_t7 = str2num(line_4(a3(3)+4 : a3(4)+3))


    T = linspace(mid_temp,max_temp,1000); 
    R = 8.314;

    Specific_heat = cp_func(R,h_t1,h_t2,h_t3,h_t4,h_t5,l_t1,l_t2,l_t3,l_t4,l_t5,mid_temp,T)
    Enthalpy = enthalpy_func(R,h_t1,h_t2,h_t3,h_t4,h_t5,h_t6,l_t1,l_t2,l_t3,l_t4,l_t5,l_t6,mid_temp,T)
    Entropy = entropy_func(R,h_t1,h_t2,h_t3,h_t4,h_t5,h_t7,l_t1,l_t2,l_t3,l_t4,l_t5,l_t7,mid_temp,T)

    %min_t_temp1 = str2num(line_1(a(1)-3:a(1)+3))
    %max_t_temp2 = str2num(line_1(a(2)-4:a(2)+3))
    %mid_t_temp3 = str2num(line_1(a(3)-4:a(3)+3))

    mkdir(species)

    % changing directory
    cd(species)
    
    % plotting
    figure(1)
    plot(T,Specific_heat)
    xlabel('Temperature')
    ylabel('Specific heat')
    title('Specific heat vs Temperature of ',sprintf('%s',species))
    saveas(1,'Specific_heat vs T.jpg')

    figure(2)
    plot(T,Entropy)
    xlabel('Temperature')
    ylabel('Entropy')
    title('Specific heat vs Entropy of ',sprintf('%s',species))
    saveas(2,'Specific_heat vs Entropy.jpg')

    figure(3)
    plot(T,Enthalpy)
    xlabel('Temperature')
    ylabel('Enthalpy')
    title('Specific heat vs Enthalpy of ',sprintf('%s',species))
    saveas(3,'Specific_heat vs Enthalpy.jpg')
    cd ..


atomic_weights = [1,12,14,16,39.94];
%     H = atomic_weights(1)
%     C = atomic_weights(2)
%     N = atomic_weights(3)
%     O = atomic_weights(4)
%     A = atomic_weights(5)
    
at_els = ['H','C','N','O','A'];
at_wt = [1,12,14,16,39.94];
mol_wt = 0;
for k = 1 : length(species)
    for j = 1 : length(at_els)
        if species(k) == at_els(j)
           mol_wt = mol_wt + at_wt(j);
           position = j;
        end
    end
    
    n = str2double(species(k))
    if n >= 1
        mol_wt = mol_wt + at_wt(position)*(n-1);
    end
    fprintf("Molecular Weight of %s = %f",species,mol_wt)
    disp(' ')
end
end

• Now we call another loop counter for 1 : 53, because the end value 53 marks the last value of the document.

Below is the function code of the Enthalpy , Entropy & Specific Heat, 

• So, we give the function command, and call the output argument which is Specific_heat / Enthalpy / Entropy which equals to the function handle by which we call the function code in the main program, inside which we provide the values for the corresponding high temperature & low temperature variables as values which we obtain by parsing the thermodat. file, which are one of the input arguments we provide to the matlab function code which is provided in the main code.
• Mid_temp is basically the local mid temperature of the species we obtain by parsing. We use this in the if-else loop to obtain the corresponding specific heat value if it is below or above the mid temp as a condition.
• R is th universal gas constant.
• T is the set of temperature values in an array from 0 to 1000 which we obtain using the linspace command.

function[Specific_heat] = cp_func(R,h_t1,h_t2,h_t3,h_t4,h_t5,l_t1,l_t2,l_t3,l_t4,l_t5,mid_temp,T)

%ccp/R = a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4

if (T <= mid_temp)
    Specific_heat = R*(l_t1 + (l_t2 .* (T)) + (l_t3 .* (T.^2)) + (l_t4 .* (T.^3)) + (l_t5 .* (T.^4))) 
else
    Specific_heat = R*(h_t1 + (h_t2 .* (T)) + (h_t3 .* (T.^2)) + (h_t4 .* (T.^3)) + (h_t5 .* (T.^4)))
end 

end

function[Entropy] = entropy_func(R,h_t1,h_t2,h_t3,h_t4,h_t5,h_t7,l_t1,l_t2,l_t3,l_t4,l_t5,l_t7,mid_temp,T)


%S/R = a1 lnT + a2 T + a3 T^2 + a4 T^3 + a5 T^4 + a7

if (T <= mid_temp)
    Entropy = R*(l_t1.*log(T) + (l_t2 .* T) + (l_t3 .* T.^2) + (l_t4 .* T.^3) + (l_t5 .* T.^4) + (l_t7))
else
    Entropy = R*(h_t1.*log(T) + (h_t2 .* T) + (h_t3 .* T.^2) + (h_t4.* T.^3) + (h_t5 .* T.^4) + (h_t7))
end 

end

function[Enthalpy] = enthalpy_func(R,h_t1,h_t2,h_t3,h_t4,h_t5,h_t7,l_t1,l_t2,l_t3,l_t4,l_t5,l_t7,mid_temp,T)


%S/R = a1 lnT + a2 T + a3 T^2 + a4 T^3 + a5 T^4 + a7

if (T <= mid_temp)
    Enthalpy = R*(l_t1 + (l_t2 .* T)./2 + (l_t3 .* T.^2)./3 + (l_t4 .* T.^3)./4 + (l_t5 .* T.^4)./5 + (l_t7./T))
else
    Enthalpy = R*(h_t1 + (h_t2 .* T)./2 + (h_t3 .* T.^2)./3 + (h_t4.* T.^3)./4 + (h_t5 .* T.^4)./5 + (h_t7./T))
end 

end