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