%% Homework 2 - Problem 2
clear all
close all
clc

% Geometric parameters (SI unit length in [m])
R = 1e-3;               % Base radius (radius of the droplet at the substrate)
h = 0.2e-3;             % Height of the droplet
H = 5 * R;              % Total height of the domain
Rc = (R^2+h^2)/(2*h);   % Radius of curvature of the droplet

% Using the steady state axysymmetric model from Matlab
% Heat and mass transfer analogy:
% Temperature -> Density
% Heat flux -> Mass flux
% D_VA -> k (as it is steady state we don't have to specify cp and rho)

% 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; H; H; 0; 0; 0; H; H];
circle = [1; 0; h-Rc; Rc];
circle_padded = [circle; zeros(length(rect) - length(circle), 1)];
gd = [rect, circle_padded];
ns = char('R1', 'C1')';
sf = 'R1-C1';
[dl, bt] = decsg(gd, sf, ns);
geometryFromEdges(model, dl);
pdegplot(dl,"EdgeLabels","on","FaceLabels","on");
% Check the edge and face label (used for the BCs)

% Generate the mesh
m1 = generateMesh(model, 'Hmax', R*0.05);
% Adjust the mesh size based the accuracy you desire and the computer you have

% Define the PDE coefficients
% https://ch.mathworks.com/help/pde/ug/pde.thermalmodel.thermalproperties.html
D_VA = 2.5e-5;                                              % Vapor diffusion coefficient in m^2/s
thermalProperties(model,"ThermalConductivity", D_VA);

% Define BCs
T_surf = 25 ;                                               % Liquid surface temperature [°C]
rho_surf = XSteam('rhoV_T',T_surf);                         % Saturation vapor density at T_surf [kg/m³] (check XSteam doc in the same folder)
T_dew = 5 ;                                                 % Far field dew point [°C]
rho_dew = XSteam('rhoV_T',T_dew);                           % Vapor density in the far field [kg/m³]

% https://ch.mathworks.com/help/pde/ug/pde.thermalmodel.thermalbc.html
thermalBC(model, 'Edge', [1,2], 'Temperature', rho_dew);    % Top and right side
thermalBC(model, 'Edge', 3, 'HeatFlux', 0);        	        % Substrate (flux normal to the BC)
thermalBC(model, 'Edge', 4, 'HeatFlux', 0);                 % Axis of symmetry (flux normal to the BC)
thermalBC(model, 'Edge', 5, 'Temperature', rho_surf);       % Droplet interface

% Solve
result = solve(model);

% Parametrize the droplet surface
theta0 = acos(R/Rc);
theta = linspace(theta0, pi/2, 30)';
r_ = Rc * cos(theta);
z_ = h - Rc + Rc * sin(theta);

% Interpolate the solution at these points
% Evaluating temperature gradient is equivalent to evaluating density
% gradient in this model

%Evaluation at the droplet surface
[rhodr, rhodz] = evaluateTemperatureGradient(result,r_,z_); 

% Plot
figure;
pdeplot(model, 'XYData', result.Temperature, 'Mesh', 'off');
xlabel('r [m]');
ylabel('z [m]');
colormap('parula');
c = colorbar;
c.Label.String = "Vapor Density Distribution (\rho_v) [kg/m³]";
axis equal
xlim([0, H]);
ylim([0,H]);