%% ME-446 Homework 6 - Problem 2
clear all
close all
clc

% 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;

%% Hsu Model - Water at 1 atm and delta_t = 0.2

% Parameters
P_l = 1.01325e5;                            % [Pa] Liquid pressure
delta_t = 0.2e-3;                           % [m] Thermal boundary layer thickness
sigma = PropsSI('I','P',P_l,'Q',0,'Water'); % [N/m] Surface tension
rho_v = PropsSI('D','P',P_l,'Q',1,'Water'); % [kg/m³] Vapor density
T_sat = PropsSI('T','P',P_l,'Q',0,'Water'); % [K] Saturation temperature
h_l = PropsSI('H','P',P_l,'Q',0,'Water');   % [J/kg] Liquid enthalpy
h_v = PropsSI('H','P',P_l,'Q',1,'Water');   % [J/kg] Vapor enthalpy
h_lv = h_v - h_l;                           % [J/kg] Heat of vaporization

% Range of wall superheat from 0 to 12°C
theta_w_min = 0;
theta_w_max = 12;
N = 5000;
theta_w = linspace(theta_w_min, theta_w_max, N)';

% Range of active cavity sizes
rc_min = (delta_t/4)*(1+sqrt(1-(12.8*sigma*T_sat)./(rho_v*h_lv*delta_t*theta_w)));
rc_max = (delta_t/4)*(1-sqrt(1-(12.8*sigma*T_sat)./(rho_v*h_lv*delta_t*theta_w)));

% Filter values based on imaginary part threshold
tolerance = 1e-9;

% Separate filters for rc_min and rc_max where imaginary part < tolerance
valid_idx_min = abs(imag(rc_min)) < tolerance;
valid_idx_max = abs(imag(rc_max)) < tolerance;

% Create separate variables for plotting
theta_w_min_real = theta_w(valid_idx_min); % Theta values for rc_min
theta_w_max_real = theta_w(valid_idx_max); % Theta values for rc_max

rc_min_real = real(rc_min(valid_idx_min)); % Real part of rc_min
rc_max_real = real(rc_max(valid_idx_max)); % Real part of rc_max

% Plot
fig1 = figure();

% Plot rc_min and rc_max only where the imaginary part is less than tolerance
semilogy(theta_w_min_real, rc_min_real*1e3, "DisplayName", "r_{c,min}")
hold on
semilogy(theta_w_max_real, rc_max_real*1e3, "DisplayName", "r_{c,max}")

xlim([theta_w_min, theta_w_max])
ylim([0.001,1])
xlabel('Wall superheat T_w - T_{sat} [°C]')
ylabel('Cavity mouth radius [mm]') 
l = legend("Location","northeast");
t = title("Saturated Water / P_{sat} = 1 atm / \delta_t = 0.2 mm");
t.FontWeight = "normal";
t.FontSize = 7;

l.FontSize = 7;
l.Box = 'off';
l.ItemTokenSize = [15,9];
l.NumColumns = 1;

% 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, 'Hsu_Base.pdf', 'ContentType', 'vector');

%% Hsu Model - Thin down the boundary layer
delta_t_arr = [0.2e-3, 0.1e-3];

% Plot
fig2 = figure();

for i = 1:length(delta_t_arr)
    delta_t = delta_t_arr(i);

    % Range of active cavity sizes
    rc_min = (delta_t/4)*(1+sqrt(1-(12.8*sigma*T_sat)./(rho_v*h_lv*delta_t*theta_w)));
    rc_max = (delta_t/4)*(1-sqrt(1-(12.8*sigma*T_sat)./(rho_v*h_lv*delta_t*theta_w)));
    
    % Filter values based on imaginary part threshold
    tolerance = 1e-9;
    
    % Separate filters for rc_min and rc_max where imaginary part < tolerance
    valid_idx_min = abs(imag(rc_min)) < tolerance;
    valid_idx_max = abs(imag(rc_max)) < tolerance;
    
    % Create separate variables for plotting
    theta_w_min_real = theta_w(valid_idx_min); % Theta values for rc_min
    theta_w_max_real = theta_w(valid_idx_max); % Theta values for rc_max
    
    rc_min_real = real(rc_min(valid_idx_min)); % Real part of rc_min
    rc_max_real = real(rc_max(valid_idx_max)); % Real part of rc_max

    % Plot rc_min and rc_max only where the imaginary part is less than tolerance
    semilogy(theta_w_min_real, rc_min_real*1e3, "DisplayName", "\delta_t = "+sprintf("%.1f mm", delta_t*1e3))
    hold on
    set(gca,'ColorOrderIndex',i)
    semilogy(theta_w_max_real, rc_max_real*1e3, "HandleVisibility", "off")

end

xlim([theta_w_min, theta_w_max])
ylim([0.001,1])
xlabel('Wall superheat T_w - T_{sat} [°C]')
ylabel('Cavity mouth radius [mm]') 
l = legend("Location","northeast");
t = title("Saturated Water / P_{sat} = 1 atm");
t.FontWeight = "normal";
t.FontSize = 7;

l.FontSize = 7;
l.Box = 'off';
l.ItemTokenSize = [15,9];
l.NumColumns = 1;

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

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

%% Hsu Model - Water / R134a
fluid_arr = ["Water", "R134a"];

% Plot
fig3 = figure();

for i = 1:length(fluid_arr)
    fluid = fluid_arr(i);

    % Parameters
    P_l = 1.01325e5;                            % [Pa] Liquid pressure
    delta_t = 0.2e-3;                           % [m] Thermal boundary layer thickness
    sigma = PropsSI('I','P',P_l,'Q',0,fluid); % [N/m] Surface tension
    rho_v = PropsSI('D','P',P_l,'Q',1,fluid); % [kg/m³] Vapor density
    T_sat = PropsSI('T','P',P_l,'Q',0,fluid); % [K] Saturation temperature
    h_l = PropsSI('H','P',P_l,'Q',0,fluid);   % [J/kg] Liquid enthalpy
    h_v = PropsSI('H','P',P_l,'Q',1,fluid);   % [J/kg] Vapor enthalpy
    h_lv = h_v - h_l;                           % [J/kg] Heat of vaporization

    % Range of active cavity sizes
    rc_min = (delta_t/4)*(1+sqrt(1-(12.8*sigma*T_sat)./(rho_v*h_lv*delta_t*theta_w)));
    rc_max = (delta_t/4)*(1-sqrt(1-(12.8*sigma*T_sat)./(rho_v*h_lv*delta_t*theta_w)));
    
    % Filter values based on imaginary part threshold
    tolerance = 1e-9;
    
    % Separate filters for rc_min and rc_max where imaginary part < tolerance
    valid_idx_min = abs(imag(rc_min)) < tolerance;
    valid_idx_max = abs(imag(rc_max)) < tolerance;
    
    % Create separate variables for plotting
    theta_w_min_real = theta_w(valid_idx_min); % Theta values for rc_min
    theta_w_max_real = theta_w(valid_idx_max); % Theta values for rc_max
    
    rc_min_real = real(rc_min(valid_idx_min)); % Real part of rc_min
    rc_max_real = real(rc_max(valid_idx_max)); % Real part of rc_max

    % Plot rc_min and rc_max only where the imaginary part is less than tolerance
    semilogy(theta_w_min_real, rc_min_real*1e3, "DisplayName", fluid)
    hold on
    set(gca,'ColorOrderIndex',i)
    semilogy(theta_w_max_real, rc_max_real*1e3, "HandleVisibility", "off")

end

xlim([theta_w_min, theta_w_max])
ylim([0.001,1])
xlabel('Wall superheat T_w - T_{sat} [°C]')
ylabel('Cavity mouth radius [mm]') 
l = legend("Location","northeast");
t = title("P_{sat} = 1 atm / \delta_t = 0.2 mm");
t.FontWeight = "normal";
t.FontSize = 7;

l.FontSize = 7;
l.Box = 'off';
l.ItemTokenSize = [15,9];
l.NumColumns = 1;

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

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

