Main [.m]

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_wall  = 1.5;    % ceramic 
-lambda_fluid = 0.6089;    % water 
+lambda_fluid = 0.6;    % water 
 lambda_air   = 0.026;  % air
-%Task 4) 
+%Task 4) heat transfer coefficient
-alpha = 10;      % heat transfer coefficient
+alpha=10;
-u_out = 18;      % ambient air temperature (°C)
+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_wall  = 1085;
-cp_fluid = 4184;
+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 
+tau = 1;              % time step in seconds 
-T_end = 400;            % total simulation time (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 
+innerWallNodes = unique([findNodes(model.Mesh,'region','Edge',8), findNodes(model.Mesh,'region','Edge',9)]);
-u = u0;
+%innerWallNodes = findNodes(model.Mesh,'region','Edge',9);
+
+u= u0;
 
 % Storage
 timeVec = (0:Nt-1)' * tau;
@@ -226,11 +249,10 @@ 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 
@@ -241,14 +263,14 @@ for k = 1:Nt
     end
 end
 
-% figure(9)
+figure(9)
-% plot(timeVec(1:k), innerWallTemp(1:k), 'LineWidth', 2)
+plot(timeVec(1:k), innerWallTemp(1:k), 'LineWidth', 2)
-% hold on
+hold on
-% yline(T_target,'r--','67°C','LineWidth',1.5)
+yline(T_target,'r--','67°C','LineWidth',1.5)
-% xlabel('Time [s]')
+xlabel('Time [s]')
-% ylabel('Average inner wall temperature [°C]')
+ylabel('Average inner wall temperature [°C]')
-% title('Heating of the inner ceramic wall')
+title('Heating of the inner ceramic wall')
-% grid on
+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;