%% Homework 9
clear all
close all
clc

% Requiring Python installed with CoolProp Library
% http://www.coolprop.org/coolprop/wrappers/Python/index.html#automatic-installation

% Define a shorter alias for the CoolProp PropsSI Python function
% See http://www.coolprop.org/coolprop/HighLevelAPI.html#propssi-function
PropsSI = @py.CoolProp.CoolProp.PropsSI;

% Define figure properties
opts.Colors     = get(groot,'defaultAxesColorOrder');
opts.width      = 8;
opts.height     = 6;
opts.fontType   = 'Arial';
opts.fontSize   = 8;

%% Problem 2 - Flow boiling heat transfer coefficient
fluid = 'R22';
P = 100e3;      % Pressure [Pa]
G = 200;        % Mass flux [kg/m²s]
q_wall = 10e3;  % Wall heat flux [W/m²]
D = 0.9e-2;     % Tube diameter [m]

% Saturation properties
rho_l = PropsSI('D','P',P,'Q',0,fluid);         % Liquid density [kg/m³]
rho_v = PropsSI('D','P',P,'Q',1,fluid);         % Vapor density [kg/m³]
mu_l = PropsSI('V','P',P,'Q',0,fluid);          % Dynamic viscosity [Pa.s]
k_l = PropsSI('L','P',P,'Q',0,fluid);           % Thermal conductivity [W/m/K]
Pr_l = PropsSI('PRANDTL','P',P,'Q',0,fluid);    % Liquid Prandtl number [-]
h_l = PropsSI('H','P',P,'Q',0,fluid);           % Liquid enthalpy [J/kg]
h_v = PropsSI('H','P',P,'Q',1,fluid);           % Vapor enthalpy [J/kg]
h_lv = h_v - h_l;                               % Heat of vaporization [J/kg]

% Non dimensional numbers
Bo = q_wall / (G * h_lv);
Re_le = (G * D) / mu_l;                 
f = (1.58 * log(Re_le) - 3.28)^(-2);
f2_Fr_le = 1;                           % For vertical tubes
Fk = 2.2;                               % for R22, from Table 12.1 in Carey

% Single phase liquid HTC
if ~(Pr_l >= 0.5 && Pr_l <= 2000)
    error("Pr_l >= 0.5 && Pr_l <= 2000 condition not respected")
end
if Re_le >= 2300 && Re_le < 1e4
    % h_le from Gninlinski (eq. 12.71 in Carey)
    h_le = (k_l / D) * (Re_le-1000)*Pr_l*(f/2) / (1 + 12.7 * (Pr_l^(2/3)-1)*(f/2)^0.5);
elseif Re_le >= 1e4 && Re_le <= 5e6
    % h_le from Petukhov and Popov (eq. 12.72 in Carey)
    h_le = (k_l / D) * Re_le*Pr_l*(f/2) / (1.07 + 12.7 * (Pr_l^(2/3)-1)*(f/2)^0.5);
end

% Vapor qualities between 20 and 60%
x_arr = linspace(0.2, 0.6, 50);
h_kandlikar_arr = [];
h_gw_arr = [];

for x = x_arr
    % Nucleate boiling regime
    h_nbd = 0.6683*(rho_l/rho_v)^0.1*x^0.16*(1-x)^0.64*f2_Fr_le*h_le+...
        1058.0*Bo^0.7*Fk*(1-x)^0.8*h_le;

    % Convective boiling regime
    h_cdb = 1.1360*(rho_l/rho_v)^0.45*x^0.72*(1-x)^0.08*f2_Fr_le*h_le+...
        667.2*Bo^0.7*Fk*(1-x)^0.8*h_le;

    % Kandlikar correlation
    h_kandlikar_arr = [h_kandlikar_arr; max(h_nbd, h_cdb)];

    % Gungor and Winterton correlation
    Re_l = G*(1-x)*D/mu_l; % Doubled checked in 1986 paper
    h_l = 0.023*(k_l/D)*Re_l^0.8*Pr_l^0.4;  % Dittus-Boetler
    h_gw = h_l*(1+3000*Bo^0.86+(x/(1-x))^0.75*(rho_l/rho_v)^0.41);
    h_gw_arr = [h_gw_arr; h_gw];
end

% Figure
fig1 = figure;
plot(x_arr, h_kandlikar_arr/1e3, DisplayName="Kandlikar");
hold on
plot(x_arr, h_gw_arr/1e3, DisplayName="Gungor & Winterton");
xlabel('Vapor Quality x [-]');
ylabel('Heat Transfer Coefficient [kW/m^2K]');
l = legend;
l.Location = "southeast";
l.Box = "off";
ylim([0, 6])

% Figure properties
fig1.Units               = 'centimeters';
fig1.Position(3)         = opts.width;
fig1.Position(4)         = opts.height;
set(fig1.Children, ...
    'FontName',     opts.fontType , ...
    'FontSize',     opts.fontSize);

% Save the plot
exportgraphics(gcf, 'HW9_P2_SOL.pdf', 'ContentType', 'vector');

%% Problem 3 - Onset of nucleate boiling
fluid = 'Nitrogen';
P = 360e3;          % Pressure [Pa]
G = 800;            % Mass flux [kg/m²s]
q_wall = 20e3;       % Wall heat flux [W/m²]
D = 0.7e-2;         % Tube diameter [m]
T_in_l = 80;        % Inlet liquid temeprature [K]

% Saturation properties
T_sat = PropsSI('T','P',P,'Q',0,fluid);         % Saturation temperature [K]
sigma = PropsSI('I','P',P,'Q',0,fluid);         % Surface tension [N/m]
cp_l = PropsSI('C','P',P,'Q',0,fluid);          % Liquid constant pressure specific heat [J/(kg.K)]
rho_v = PropsSI('D','P',P,'Q',1,fluid);         % Vapor density [kg/m³]
mu_l = PropsSI('V','P',P,'Q',0,fluid);          % Dynamic viscosity [Pa.s]
Pr_l = PropsSI('PRANDTL','P',P,'Q',0,fluid);    % Liquid Prandtl number [-]
k_l = PropsSI('L','P',P,'Q',0,fluid);           % Thermal conductivity [W/m/K]
h_l = PropsSI('H','P',P,'Q',0,fluid);           % Liquid enthalpy [J/kg]
h_v = PropsSI('H','P',P,'Q',1,fluid);           % Vapor enthalpy [J/kg]
h_lv = h_v - h_l;                               % Heat of vaporization [J/kg]

% Single phase liquid HTC
Re_le = (G*D)/mu_l;
h_le = 0.023*(k_l/D)*Re_le^0.8*Pr_l^0.4;  % Dittus-Boetler

% Needed wall superheat for the ONB
delta_T_w = sqrt((q_wall*8*sigma*T_sat)/(k_l*h_lv*rho_v));

% Subcooling
delta_T_sub = T_sat - T_in_l;


% Position of the ONB
z_onb = (G*cp_l*D)/(4*h_le)*((h_le/q_wall)*(delta_T_w+delta_T_sub)-1);
fprintf("z_ONB = %3f m\n", z_onb)