From 9cb746341c741691110d81a98170b9570f093ee5 Mon Sep 17 00:00:00 2001 From: "lisa.pizzo" Date: Sun, 25 Jan 2026 11:49:32 +0100 Subject: [PATCH] Main file [.m] --- Project/Main.m | 298 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 298 insertions(+) create mode 100644 Project/Main.m diff --git a/Project/Main.m b/Project/Main.m new file mode 100644 index 0000000..23585ed --- /dev/null +++ b/Project/Main.m @@ -0,0 +1,298 @@ +% Axisymmetric mug +% r–z plane, rotation around r = 0 +clear; clc; close all; + +%% General values that we use in the entire script +maxIter = 5000; +tol = 1e-6; +%Task 3) Thermal conductivity +lambda_wall = 2.0; % ceramic RANDOM NUMBER +lambda_fluid = 0.6; % water RANDOM NUMBER +lambda_air = 0.025; % air RANDOM NUMBER +%Task 4) +alpha = 10; % heat transfer coefficient +u_out = 18; % ambient air temperature (°C) +%Task 6) Heat capacities (example values) +c_wall = 900; % ceramic RANDOM NUMBER +c_fluid = 4180; % water RANDOM NUMBER +c_air = 1005; % air RANDOM NUMBER + +% This numbers don't simulate real life because in real life we have +% volumetric heat capacities. this means we are having a reduced factore of +% approx 10^4. + +%% Task 1: Mesh definition +% Create PDE model +model = createpde(); + +% Points (meters) +A = [0, 0]; +B = [0.055, 0]; +C = [0.083, 0.105]; +H = [0.078, 0.105]; +F = [0.050, 0.005]; +E = [0, 0.005]; +G = [0.067, 0.066]; +I = [0, 0.066]; +D = [0, 0.105]; + +% Geometry matrix (edges) +% Axis (split by material) +g1 = [2; A(1); E(1); A(2); E(2); 1; 0]; % ceramic +g2 = [2; E(1); I(1); E(2); I(2); 2; 0]; % fluid +g3 = [2; I(1); D(1); I(2); D(2); 3; 0]; % air +% Outer ceramic +g4 = [2; A(1); B(1); A(2); B(2); 1; 0]; +g5 = [2; B(1); C(1); B(2); C(2); 1; 0]; +% Top rim: C -> H is ceramic-air (ONLY) +g6 = [2; C(1); H(1); C(2); H(2); 1; 3]; +% Inner ceramic wall (H -> F) ceramic-air +g7 = [2; H(1); F(1); H(2); F(2); 1; 3]; +% Inner ceramic bottom (F -> E) ceramic-fluid +g8 = [2; F(1); E(1); F(2); E(2); 1; 2]; +% Fluid surface (F -> G) fluid-air +g9 = [2; F(1); G(1); F(2); G(2); 2; 3]; +% Fluid surface (G -> I) fluid-air +g10 = [2; G(1); I(1); G(2); I(2); 2; 3]; +% Air top boundary: D -> H (air-outside) +g11 = [2; D(1); H(1); D(2); H(2); 3; 0]; + +% Assemble geometry +g = [g1 g2 g3 g4 g5 g6 g7 g8 g9 g10 g11]; + +% Load geometry into PDE model +geometryFromEdges(model, g); + +% Plot geometry with labels +% figure(1); +% pdegplot(model, 'EdgeLabels','on', 'FaceLabels','on'); +% axis equal; +% title('Geometry with edge and face labels'); + +% Generate mesh +mesh = generateMesh(model, 'Hmax', 0.006, 'GeometricOrder','linear'); + +% Plot mesh +% figure(2); +% pdemesh(model); +% axis equal; +% title('Generated mesh'); + +%% Task 2: Jacobi solver with constant lambda +nodes = mesh.Nodes; % coordinates of all mesh nodes +elements = mesh.Elements; % which nodes make up each triangle element + +Nnodes = size(nodes,2); +Nelems = size(elements,2); + +% Define material properties (for simplicity, lambda = 1 everywhere) +lambda = ones(Nelems,1); % thermal conductivity + +% Initialize global stiffness matrix and RHS +K = sparse(Nnodes, Nnodes); +F = zeros(Nnodes,1); + +% Assemble K and F +for e = 1:Nelems %Loop over each triangle element + vert = elements(:,e); %nodes of element + x = nodes(1,vert); + y = nodes(2,vert); + + Ae = polyarea(x,y); % Compute area of the triangle + + % Linear triangle gradients + b = [y(2)-y(3); y(3)-y(1); y(1)-y(2)]; % derivative with respect to x + c = [x(3)-x(2); x(1)-x(3); x(2)-x(1)]; % derivative with respect to y + + % Element stiffness matrix + Ke = (lambda(e)/(4*Ae)) * (b*b.' + c*c.'); + + % Assemble + K(vert,vert) = K(vert,vert) + Ke; + + % Element load vector (f=0) + F(vert) = F(vert) + zeros(3,1); +end + +% Find boundary nodes (couldn't find a better way) +edgesAll = [elements([1 2],:), elements([2 3],:), elements([3 1],:)]; % all edges -> 2x(3*Nelems) awway since each triangle has 3 edges +edgesSorted = sort(edgesAll,1); % sort nodes of each edge, ensure that [i,j] and [j,i] are the same +[~,~,ic] = unique(edgesSorted','rows'); % identifies unique edges and assignes indices ic +counts = accumarray(ic,1); % counts how many times each edge appears in the mesh +boundaryEdges = find(counts==1); %edges belonging to only 1 element: hence on the boundary +boundaryNodes = unique(edgesSorted(:,boundaryEdges)); % nodes belonging to these boundary edges + +% Jacobi iteration +D = diag(K); +R = K - diag(D); +u = zeros(Nnodes,1); + +for iter = 1:maxIter + u_new = (F - R*u) ./ D; + u_new(boundaryNodes) = 0; % Enforce Dirichlet BC + if norm(u_new - u, inf) < tol % Convergence check + fprintf('Task 2) Jacobi converged in %d iterations.\n', iter); + break; + end + u = u_new; +end + +% % Plot solution +% figure(3) +% pdeplot(model, 'XYData', u, 'Mesh','on'); +% axis equal; +% title('Stationary Dirichlet solution via Jacobi'); +% colorbar; + +%% Task 3: Laplace with multiple lambdas +[K, F] = CalculateLaplace_mult(model, lambda_wall, lambda_fluid, lambda_air); + +D = diag(K); +R = K - diag(D); +u = zeros(Nnodes,1); +for iter = 1:maxIter + u_new = (F - R*u) ./ D; + u_new(boundaryNodes) = 0; % Dirichlet BC + + if norm(u_new - u, inf) < tol + fprintf('Task 3) Jacobi converged in %d iterations.\n', iter); + break; + end + u = u_new; +end + +% figure(4) +% pdeplot(model, 'XYData', u, 'Mesh','on'); +% axis equal +% title('Task 3: Stationary solution with multiple conductivities'); +% colorbar + +%% Task 4: Robin boundary condition +[K, F] = ApplyRobinBC_mult(model, K, F, alpha, u_out); +D = diag(K); +R = K - diag(D); +u = zeros(Nnodes,1); + +for iter = 1:maxIter + u_new = (F - R*u) ./ D; + + if norm(u_new - u, inf) < tol + fprintf('Task 4) Jacobi (Robin BC) converged in %d iterations.\n', iter); + break; + end + u = u_new; +end + +figure(5) +pdeplot(model, 'XYData', u, 'Mesh','on'); +axis equal +title('Task 4: Stationary solution with Robin BC'); +colorbar + +%% 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); + +D = diag(K); +R = K - diag(D); +u = zeros(Nnodes,1); + +for iter = 1:maxIter + u_new = (F - R*u) ./ D; + if norm(u_new - u, inf) < tol + fprintf('Task 5) Jacobi (axisymmetric Robin) converged in %d iterations.\n', iter); + break; + end + u = u_new; +end + +% Plot solution +% figure(6) +% pdeplot(model, 'XYData', u, 'Mesh','on'); +% axis equal +% title('Task 5: Axisymmetric stationary solution with Robin BC'); +% colorbar + +%To see it in 3D paste here the code in "AdditionalPlotCodes.txt". + +%% Task 6: Axisymmetric mass matrix +M = sparse(Nnodes, Nnodes); +M = AddMass_mult_rot(model, M, c_wall, c_fluid, c_air); + +%% Task 7: Initial solution +u0 = Init_Solution_mult(model, 18, 80, 18); + +% figure(7) +% pdeplot(model, 'XYData', u0, 'Mesh','on'); +% axis equal +% title('Initial temperature distribution'); +% colorbar + +%% Task 8: Time-dependent simulation (explicit scheme) +tau = 0.5; % time step in seconds +T_end = 400; % total simulation time (seconds) +Nt = ceil(T_end/tau); % number of time steps + +A = (1/tau)*M+K; % Left-hand side matrix +[L,U,P,Q] = lu(A); % A is constant -> factorize it once: PAQ=LU + +% Initialize solution +u = u0; + +for k = 1:Nt + b = (1/tau)*M*u + F; % F is the load vector, F=0 + + % Solve for next time step + u_next = Q*(U\(L\(P*b))); % efficient solution using LU factors + + u = u_next; % Update + + if mod(k,20) == 0 + % figure(8) + % pdeplot(model, 'XYData', u, 'Mesh','on'); + % axis equal + % title(['Temperature at t = ', num2str(k*tau), ' s']); + % colorbar + % drawnow + end +end + +%To see the 9 snapshots paste here the codes in "AdditionalPlotCodes.txt" + +%% Task 9 (i): Heating time using inner ceramic wall temperature +T_target = 67; % [°C] + +innerWallNodes = findNodes(model.Mesh,'region','Edge',8); % Edge 8 = ceramic–fluid interface +u = u0; + +% Storage +timeVec = (0:Nt-1)' * tau; +innerWallTemp = zeros(Nt,1); +Twarm = NaN; + +for k = 1:Nt + b = (1/tau)*M*u + F; + + u = Q*(U\(L\(P*b))); + + % Average inner wall temperature + innerWallTemp(k) = mean(u(innerWallNodes)); + + % Check heating criterion + 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); + break + end +end + +%% --- Plot inner wall temperature evolution +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