%% ME-446 Homeowrk 10 
clear all
close all
clc

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

%% Problem 1 - Thermal conduction in a droplet

% Geometric parameters (SI unit length in [m])
d = 5e-3;               % Droplet diameter (radius of the droplet at the substrate)
r = d/2;

% Using the steady state axysymmetric model from Matlab

% Create a new PDE model
% https://ch.mathworks.com/help/pde/ug/createpde.html
model = createpde('thermal', 'steadystate-axisymmetric');

% Generate geometry
% https://ch.mathworks.com/help/pde/ug/decsg.html
% https://ch.mathworks.com/help/pde/ug/pde.pdemodel.geometryfromedges.html
rect = [3; 4; 0; r; r; 0; 0; 0; r; r];
circle = [1; 0; 0; r];
circle_padded = [circle; zeros(length(rect) - length(circle), 1)];
gd = [rect, circle_padded];
ns = char('R1', 'C1')';
sf = 'C1*R1';
[dl, bt] = decsg(gd, sf, ns);
geometryFromEdges(model, dl);

% Plot
fig1 = figure(1);
pdegplot(dl,"EdgeLabels","on","FaceLabels","on");

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

% Generate the mesh
N = 100;
mesh = generateMesh(model, 'Hmax', r/N);

% Savings the Nodes on Edge 1 and 3
Ne1 = findNodes(mesh,"region","Edge",1);
Ne3 = findNodes(mesh,"region","Edge",3);

% Define the PDE coefficients
% https://ch.mathworks.com/help/pde/ug/pde.thermalmodel.thermalproperties.html
k = 0.6;     % Thermal conductivity [W/m.K]
thermalProperties(model,"ThermalConductivity", k);

% Define BCs
Ta = 50 ;                         % Interface temperature [°C]
h = 10e3;                           % Heat transfer coefficient at the interface [W/(m².K)]
Ts_arr = linspace(10, 40, 15)' ;    % Substrate temperature [°C]

thermalBC(model, 'Edge', 2, 'HeatFlux', 0);        	                                    % Axis of symmetry (zero flux normal to the BC)
thermalBC(model, 'Edge', 3, 'ConvectionCoefficient', h, 'AmbientTemperature', Ta);      % Interface

% Declaring zeros arrays to store the heat flow rate across the droplet qd 
% and interface temperature Ti
qd = zeros(length(Ts_arr),1);    
Ti = zeros(length(Ts_arr),1);

% Loop over the substrate temperature array
for i = 1:length(Ts_arr)

    % Apply the temperature BC at the substrate
    t_s = Ts_arr(i);
    thermalBC(model, 'Edge', 1, 'Temperature', t_s);
    
    % Solve
    result = solve(model);
    
    % Evaluation of the heat rate at the droplet-substrate interface
    dTdz_s = result.ZGradients(Ne1);
    r_s = result.Mesh.Nodes(1,Ne1)';    % First row of result.Mesh.Nodes gives r-coordinates of the nodes
   
    % Heat flux along the z direction
    qz = -k*dTdz_s;

    % Integrate over the radial component from 0 to r and axial component from
    % 0 to 2π. Taking only the z component of the heat flux as we look for the
    % flow rate going out.
    q(i) = abs(2*pi*trapz(r_s, r_s.*qz)); % [W]

    % Note that we obtain a similar result using the built-in
    % evaluateHeatRate function
    q(i) = evaluateHeatRate(result, Edge=1);
    
    % Average temperature at the interface Ti (needed to calculate DeltaT_cond)
    Ti(i) = mean(result.Temperature(Ne3));

end

DeltaT_cond = Ti - Ts_arr;

% Plot
fig2 = figure(2);
plot(DeltaT_cond, q, "DisplayName", "Numerical")
hold on

% Eq. 9.28 from Carey
q_carey = DeltaT_cond*2*pi*k*r;
h = plot(DeltaT_cond, q_carey, "DisplayName", "Carey Eq. 9.28");
h.LineStyle = "--";

xlabel('\Delta T_{cond} [°C]');
ylabel('Heat flow rate [W]');
leg = legend("Location", "northwest");

leg.FontSize = 7;
leg.Box = 'off';
leg.ItemTokenSize = [15,9];
leg.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, 'HW10_P1_FIG_SOL.pdf', 'ContentType', 'vector');

%% Problem 2 - Dropwise condensation correlation
clear all
close all
clc

% Import the CoolProp module
CP = py.importlib.import_module('CoolProp.CoolProp');
f = 'Water';

% Parameters
t = table();
t.p = linspace(50, 9460, 200)';   % [kPa] Range of pressure
delta_t_sub = 5;                 % [°C] Subcooling

% Preallocate columns
t.h_dc = zeros(height(t), 1);
t.k_l = zeros(height(t), 1);
t.sigma = zeros(height(t), 1);
t.mu_l = zeros(height(t), 1);

% Loop over the pressures and calculate properties
for i = 1:height(t)
    p = t.p(i) * 1000;  % Convert kPa to Pa

    % Get the fluid properties at saturation
    t_v = CP.PropsSI('T', 'P', p, 'Q', 1, f);            % Saturation temperature [K]
    t.k_l(i) = CP.PropsSI('L', 'P', p, 'Q', 0, f);       % Liquid thermal conductivity [W/m.K]
    t.sigma(i) = CP.PropsSI('I', 'P', p, 'Q', 0, f) * 1e3; % Surface tension [mN/m]
    t.mu_l(i) = CP.PropsSI('V', 'P', p, 'Q', 0, f) * 1e6; % Dynamic viscosity [µPa.s]

    % Rose et al. Correlation (t_v in [°C], and h_dc in [kW/m²K])
    t.h_dc(i) = (t_v - 273.15)^0.8 * (5 + 0.3 * delta_t_sub);
end

% Plotting with subplot
figure;

% Heat transfer coefficient
subplot(4, 1, 1);
plot(t.p, t.h_dc);
ylabel('h_{dc} [kW/m²K]');
grid on;
set(gca, 'XTickLabel', []); % Remove x-axis tick labels
set(gca, 'XTick', []);      % Remove x-axis ticks

% Liquid thermal conductivity
subplot(4, 1, 2);
plot(t.p, t.k_l);
ylabel('k_l [W/(m.K)]');
grid on;
set(gca, 'XTickLabel', []); % Remove x-axis tick labels
set(gca, 'XTick', []);      % Remove x-axis ticks

% Surface tension
subplot(4, 1, 3);
plot(t.p, t.sigma); 
ylabel('\sigma [N/m]');
grid on;
set(gca, 'XTickLabel', []); % Remove x-axis tick labels
set(gca, 'XTick', []);      % Remove x-axis ticks

% Dynamic viscosity
subplot(4, 1, 4);
plot(t.p, t.mu_l);
xlabel('Pressure [kPa]');
ylabel('\mu_l [µPa.s]');
grid on;

%Tight layout adjustment
set(gcf, 'Position', [100, 100, 600, 800]); % Adjust figure size

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