%% ME-446 Homeowrk 9 - Problem 1
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;

% Geometrical Parameters
l_arr_ = [linspace(4, 60, 200)',...
         linspace(8, 60, 200)',...
         linspace(11, 60, 200)',...
         linspace(15, 60, 200)']* 1e-6;         % Micropillars pitch [m]
d_arr = [2.5, 5, 7.5, 10]' * 1e-6;              % Micropillars diameter [m]
h = 25e-6;                                      % Micropillars height [m]
W = 500e-6;                                     % Microchannel width [m]
Lw = W;                                         % Longuest wicking distance [m]

% Fluid properties
mu = 281.58e-6;                                 % Liquid viscosity @100°C [Pa.s]
gamma = 58.912e-3;                              % Surface tension @100°C [N/m]

% Surface interaction
theta_r = deg2rad(15);

% Plot
fig1 = figure(1);

% Loop over all the geometries
for i = 1:length(d_arr)
    d = d_arr(i);
    l_arr = l_arr_(:,i);
    u_ave = zeros(length(l_arr),1);

    for j = 1:length(l_arr)
        l = l_arr(j);

        % Calculate porosity and permeability
        c = (pi*d^2)/(4*l^2);                                       % Solid fraction [-]
        epsilon = 1 - c;                                            % Porosity [-]
        K = l^2*(log(c^(-1/2))-0.738+c-0.887*c^2+2.038*c^3)/(4*pi); % Permeability [m²]
        alpha = 1/sqrt(K);
        
        
        % Capilary pressure
        DeltaP_cap = (gamma*pi*d*cos(theta_r))/(l^2-0.25*pi*d^2);
        dPdx = -DeltaP_cap / Lw;
        
        % Constants A and B
        A = dPdx*exp(-alpha*h*sqrt(epsilon))/(alpha^2*mu*(exp(alpha*h*sqrt(epsilon))+exp(-alpha*h*sqrt(epsilon))));
        B = dPdx*exp(alpha*h*sqrt(epsilon))/(alpha^2*mu*(exp(alpha*h*sqrt(epsilon))+exp(-alpha*h*sqrt(epsilon))));

        % Analytical solution of the integral
        a = alpha*sqrt(epsilon);
        C = (1/(alpha^2*mu))*dPdx;
        u_ave(j) = (1/h)*((A/a)*(exp(a*h)-1)+(B/a)*(1-exp(-a*h))-C*h);

    end
    
    plot(l_arr*1e6, u_ave, "DisplayName", sprintf("d = %.1f um", d*1e6))
    hold on
end

xlabel('Pillar Pitch l [um]')
ylabel('Liquid Velocity u_{ave} [m/s]') 
leg = legend("Location","northwest");

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

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

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