%% Homework 6
clear all
close all
clc

% Problem 1 - Rohsenow's Correlation

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


% Using Coolprop's high-level interface to get the fluid properties
fluid = 'Water';
P_l = 1.01325e5;                                    % Liquid Pressure [Pa]
T_sat = PropsSI('T','P',P_l,'Q',0,fluid);           % Saturation temperature [K]
mu_l = PropsSI('V','P',P_l,'Q',0,fluid);            % Dynamic viscosity [Pa.s]
rho_l = PropsSI('D','P',P_l,'Q',0,fluid);           % Liquid density [kg/m³]
rho_v = PropsSI('D','P',P_l,'Q',1,fluid);           % Vapor density [kg/m³]
sigma = PropsSI('I','P',P_l,'Q',0,fluid);           % Surface tension [N/m]
h_l = PropsSI('H','P',P_l,'Q',0,fluid);             % Liquid enthalpy [J/kg]
h_v = PropsSI('H','P',P_l,'Q',1,fluid);             % Vapor enthalpy [J/kg]
h_lv = h_v - h_l;                                   % Heat of vaporization [J/kg]
cp_l = PropsSI('C','P',P_l,'Q',0,fluid);            % Liquid constant pressure specific heat [J/(kg.K)]
pr_l = PropsSI('PRANDTL','P',P_l,'Q',0,fluid);      % Liquid Prandtl number [-]
g = 9.81;                                           % Gravitational acceleration [m/s²]

% Rohsenow constants for water on mechanically polished stainless steel
r = 0.33;
s = 1;
c_sf = 0.0132;

% Wall superheat range
T_w = linspace(110, 130, 100)' + 273.15;

% Heat flux / Rohsenow correlation
q = mu_l*h_lv/(sigma/(g*(rho_l-rho_v)))^(0.5) * (1/c_sf)^(1/r) * ...
    pr_l^(-s/r) * (cp_l*(T_w-T_sat)/h_lv).^(1/r);

% Boiling curve
fig1 = figure();
plot(T_w-T_sat, q/1e3)
xlim([10, 30])
xlabel('T_{w}-T_{sat}(P_l) [K]')
ylabel('q" [kW/m²]') 

t = title("Rohsenow correlation for Water at P_{l} = 1 atm");
t.FontWeight = "normal";
t.FontSize = 7;

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