From 571d31823543f0bd87898d453f8591f597c5617f Mon Sep 17 00:00:00 2001 From: "lisa.pizzo" Date: Wed, 28 Jan 2026 13:13:13 +0100 Subject: [PATCH] Main [.m] --- Project/Main.m | 78 ++++++++++++++++++++++++++++++++------------------ 1 file changed, 50 insertions(+), 28 deletions(-) diff --git a/Project/Main.m b/Project/Main.m index fa32fe4..ca0461e 100644 --- a/Project/Main.m +++ b/Project/Main.m @@ -4,20 +4,20 @@ clear; clc; close all; %% General values that we use in the entire script %Task 3) Thermal conductivity -lambda_wall = 30; % ceramic -lambda_fluid = 0.6089; % water +lambda_wall = 1.5; % ceramic +lambda_fluid = 0.6; % water lambda_air = 0.026; % air -%Task 4) -alpha = 10; % heat transfer coefficient -u_out = 18; % ambient air temperature (°C) +%Task 4) heat transfer coefficient +alpha=10; +u_out = 18; % ambient air temperature (°C) %Task 6) Volumetric heat capacities: rho * c_p % Densities [kg/m^3] rho_wall = 2400; % ceramic rho_fluid = 1000; % water rho_air = 1.2; % air % Specific heats [J/(kg*K)] -cp_wall = 900; -cp_fluid = 4184; +cp_wall = 1085; +cp_fluid = 4186; cp_air = 1005; % Volumetric heat capacities [J/(m^3*K)] c_wall = rho_wall * cp_wall; @@ -44,7 +44,7 @@ g5 = [2; B(1); C(1); B(2); C(2); 1; 0]; % Outer ceramic g6 = [2; C(1); H(1); C(2); H(2); 1; 3]; % Top rim: (C-H) ceramic-air g7 = [2; H(1); F(1); H(2); F(2); 1; 3]; % Inner ceramic wall (H-F) ceramic-air g8 = [2; F(1); E(1); F(2); E(2); 1; 2]; % Inner ceramic bottom (F-E) ceramic-fluid -g9 = [2; F(1); G(1); F(2); G(2); 2; 3]; % Fluid surface (F-G) fluid-air +g9 = [2; F(1); G(1); F(2); G(2); 2; 3]; % Fluid surface (F-G) fluid - ceramic g10 = [2; G(1); I(1); G(2); I(2); 2; 3]; % Fluid surface (G-I) fluid-air g11 = [2; D(1); H(1); D(2); H(2); 3; 0]; % Air top boundary: (D-H) air-outside @@ -55,7 +55,19 @@ geometryFromEdges(model, g); % figure(1); % pdegplot(model, 'EdgeLabels','on', 'FaceLabels','on'); % axis equal; -% title('Geometry with edge and face labels'); +% hold on; +% +% % Plot points +% Pts = [A; B; C; D; E; F; G; H; I]; +% plot(Pts(:,1), Pts(:,2), 'ro', 'MarkerFaceColor','r'); +% % Label points +% labels = {'A','B','C','D','E','F','G','H','I'}; +% for k = 1:length(labels) +% text(Pts(k,1), Pts(k,2), [' ' labels{k}], ... +% 'FontSize',10, 'Color','r', 'FontWeight','bold'); +% end +% title('Geometry with edge, face, and point labels'); +% hold off; % Generate mesh, linear and 3 nodes per element mesh = generateMesh(model, 'Hmax', 0.002, 'GeometricOrder','linear'); @@ -159,7 +171,7 @@ u = K \ F; %% Task 5: Axisymmetric Laplace + Robin BC [K, F] = CalculateLaplace_mult_rot(model, lambda_wall, lambda_fluid, lambda_air); -[K, F] = ApplyRobinBC_mult_rot(model, K, F, alpha, u_out); +[K,F] = ApplyRobinBC_mult_rot(model, K, F, alpha, u_out); % Direct solve u = K \ F; @@ -178,6 +190,7 @@ M = AddMass_mult_rot(model, M, c_wall, c_fluid, c_air); %% Task 7: Initial solution u0 = Init_Solution_mult(model, 18, 80, 18); +% TRY WITH HOTTER LIQUID % figure(7) % pdeplot(model, 'XYData', u0, 'Mesh','on'); @@ -186,8 +199,8 @@ u0 = Init_Solution_mult(model, 18, 80, 18); % colorbar %% Task 8: Time-dependent simulation (explicit scheme) -tau = 0.5; % time step in seconds -T_end = 400; % total simulation time (seconds) +tau = 1; % time step in seconds +T_end = 1000; % total simulation time (seconds) Nt = ceil(T_end/tau); % number of time steps A = (1/tau)*M+K; % Left-hand side matrix @@ -209,15 +222,25 @@ for k = 1:Nt end %To see the 9 snapshots paste here the codes in "AdditionalPlotCodes.txt" +M = sparse(Nnodes, Nnodes); +M = AddMass_mult_rot(model, M, c_wall, c_fluid, c_air); %% Task 9 (i): Heating time using inner ceramic wall temperature T_target = 67; % [°C] -[K, F] = CalculateLaplace_mult_rot(model, lambda_wall, lambda_fluid, lambda_air); -[K, F] = ApplyRobinBC_mult_rot(model, K, F, alpha, u_out); + +tau = 0.1; % time step in seconds +T_end = 1000; % total simulation time (seconds) +Nt = ceil(T_end/tau); % number of time steps + +[K,F] = CalculateLaplace_mult_rot(model, lambda_wall, lambda_fluid, lambda_air); +[K,F] = ApplyRobinBC_mult_rot(model, K, F, alpha, u_out); + A = (1/tau)*M+K; % Left-hand side matrix -innerWallNodes = findNodes(model.Mesh,'region','Edge',8); % Edge 8 = ceramic–fluid interface -u = u0; +innerWallNodes = unique([findNodes(model.Mesh,'region','Edge',8), findNodes(model.Mesh,'region','Edge',9)]); +%innerWallNodes = findNodes(model.Mesh,'region','Edge',9); + +u= u0; % Storage timeVec = (0:Nt-1)' * tau; @@ -226,14 +249,13 @@ Twarm = NaN; for k = 1:Nt b = (1/tau)*M*u + F; - u = A\b; - % Average inner wall temperature innerWallTemp(k) = mean(u(innerWallNodes)); + %innerWallTemp(k) = max(u(innerWallNodes)); % Check heating criterion - if innerWallTemp(k) >= T_target + if innerWallTemp(k) >= T_target Twarm = k * tau; fprintf('Task 9 (i): Inner wall reaches %.1f°C at T = %.1f s\n', ... T_target, Twarm); @@ -241,14 +263,14 @@ for k = 1:Nt end end -% figure(9) -% plot(timeVec(1:k), innerWallTemp(1:k), 'LineWidth', 2) -% hold on -% yline(T_target,'r--','67°C','LineWidth',1.5) -% xlabel('Time [s]') -% ylabel('Average inner wall temperature [°C]') -% title('Heating of the inner ceramic wall') -% grid on +figure(9) +plot(timeVec(1:k), innerWallTemp(1:k), 'LineWidth', 2) +hold on +yline(T_target,'r--','67°C','LineWidth',1.5) +xlabel('Time [s]') +ylabel('Average inner wall temperature [°C]') +title('Heating of the inner ceramic wall') +grid on %% CHECK: Insulated mug – transient redistribution [K_neu, F_neu] = CalculateLaplace_mult_rot(model,lambda_wall,lambda_fluid,lambda_air); @@ -258,7 +280,7 @@ M_neu = AddMass_mult_rot(model,M_neu,c_wall,c_fluid,c_air); % Time stepping parameters tau = 0.5; -T_end = 400; +T_end = 100; Nt = ceil(T_end / tau); A_neu = (1/tau) * M_neu + K_neu;