From 70dc0861cb1c8680fe47d4f50cb79372640c9174 Mon Sep 17 00:00:00 2001 From: "lisa.pizzo" Date: Mon, 26 Jan 2026 18:36:01 +0100 Subject: [PATCH] Delete Project/Main.m --- Project/Main.m | 298 ------------------------------------------------- 1 file changed, 298 deletions(-) delete mode 100644 Project/Main.m diff --git a/Project/Main.m b/Project/Main.m deleted file mode 100644 index 23585ed..0000000 --- a/Project/Main.m +++ /dev/null @@ -1,298 +0,0 @@ -% 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