Main [.m]
This commit is contained in:
parent
f61f29a2de
commit
287dbfb0f6
1 changed files with 281 additions and 0 deletions
281
Project/Main.m
Normal file
281
Project/Main.m
Normal file
|
|
@ -0,0 +1,281 @@
|
||||||
|
% Axisymmetric mug
|
||||||
|
% r–z plane, rotation around r = 0
|
||||||
|
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_air = 0.026; % air
|
||||||
|
%Task 4)
|
||||||
|
alpha = 10; % heat transfer coefficient
|
||||||
|
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_air = 1005;
|
||||||
|
% Volumetric heat capacities [J/(m^3*K)]
|
||||||
|
c_wall = rho_wall * cp_wall;
|
||||||
|
c_fluid = rho_fluid * cp_fluid;
|
||||||
|
c_air = rho_air * cp_air;
|
||||||
|
|
||||||
|
%% 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)
|
||||||
|
g1 = [2; A(1); E(1); A(2); E(2); 1; 0]; % Axis - ceramic
|
||||||
|
g2 = [2; E(1); I(1); E(2); I(2); 2; 0]; % Axis - fluid
|
||||||
|
g3 = [2; I(1); D(1); I(2); D(2); 3; 0]; % Axis - air
|
||||||
|
g4 = [2; A(1); B(1); A(2); B(2); 1; 0]; % Outer ceramic
|
||||||
|
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
|
||||||
|
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
|
||||||
|
|
||||||
|
% Assemble geometry
|
||||||
|
g = [g1 g2 g3 g4 g5 g6 g7 g8 g9 g10 g11];
|
||||||
|
geometryFromEdges(model, g);
|
||||||
|
|
||||||
|
% figure(1);
|
||||||
|
% pdegplot(model, 'EdgeLabels','on', 'FaceLabels','on');
|
||||||
|
% axis equal;
|
||||||
|
% title('Geometry with edge and face labels');
|
||||||
|
|
||||||
|
% Generate mesh, linear and 3 nodes per element
|
||||||
|
mesh = generateMesh(model, 'Hmax', 0.002, 'GeometricOrder','linear');
|
||||||
|
|
||||||
|
% figure(2);
|
||||||
|
% pdemesh(model);
|
||||||
|
% axis equal;
|
||||||
|
% title('Generated mesh');
|
||||||
|
|
||||||
|
%% Task 2: Direct 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
|
||||||
|
|
||||||
|
% Direct solver
|
||||||
|
% Enforce Dirichlet BC strongly
|
||||||
|
K(boundaryNodes,:) = 0;
|
||||||
|
K(:,boundaryNodes) = 0;
|
||||||
|
K(boundaryNodes,boundaryNodes) = speye(length(boundaryNodes));
|
||||||
|
F(boundaryNodes) = 0;
|
||||||
|
|
||||||
|
% Direct solve
|
||||||
|
u = K \ F;
|
||||||
|
|
||||||
|
% figure(3)
|
||||||
|
% pdeplot(model, 'XYData', u, 'Mesh','on');
|
||||||
|
% axis equal;
|
||||||
|
% title('Stationary Dirichlet solution');
|
||||||
|
% colorbar;
|
||||||
|
|
||||||
|
%% Task 3: Laplace with multiple lambdas
|
||||||
|
[K, F] = CalculateLaplace_mult(model, lambda_wall, lambda_fluid, lambda_air);
|
||||||
|
|
||||||
|
% Enforce Dirichlet BC strongly
|
||||||
|
K(boundaryNodes,:) = 0;
|
||||||
|
K(:,boundaryNodes) = 0;
|
||||||
|
K(boundaryNodes,boundaryNodes) = speye(length(boundaryNodes));
|
||||||
|
F(boundaryNodes) = 0;
|
||||||
|
|
||||||
|
% Direct solve
|
||||||
|
u = K \ F;
|
||||||
|
|
||||||
|
% % Plot solution
|
||||||
|
% 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] = CalculateLaplace_mult(model, lambda_wall, lambda_fluid, lambda_air);
|
||||||
|
[K, F] = ApplyRobinBC_mult(model, K, F, alpha, u_out);
|
||||||
|
|
||||||
|
% Direct solve
|
||||||
|
u = K \ F;
|
||||||
|
|
||||||
|
% 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);
|
||||||
|
|
||||||
|
% Direct solve
|
||||||
|
u = K \ F;
|
||||||
|
|
||||||
|
% 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
|
||||||
|
u = u0;
|
||||||
|
|
||||||
|
for k = 1:Nt
|
||||||
|
b = (1/tau)*M*u + F; % F is the load vector, F=0
|
||||||
|
u_next = A\b;
|
||||||
|
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]
|
||||||
|
[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;
|
||||||
|
|
||||||
|
% Storage
|
||||||
|
timeVec = (0:Nt-1)' * tau;
|
||||||
|
innerWallTemp = zeros(Nt,1);
|
||||||
|
Twarm = NaN;
|
||||||
|
|
||||||
|
for k = 1:Nt
|
||||||
|
b = (1/tau)*M*u + F;
|
||||||
|
|
||||||
|
u = A\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
|
||||||
|
|
||||||
|
% 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);
|
||||||
|
|
||||||
|
M_neu = sparse(Nnodes, Nnodes);
|
||||||
|
M_neu = AddMass_mult_rot(model,M_neu,c_wall,c_fluid,c_air);
|
||||||
|
|
||||||
|
% Time stepping parameters
|
||||||
|
tau = 0.5;
|
||||||
|
T_end = 400;
|
||||||
|
Nt = ceil(T_end / tau);
|
||||||
|
|
||||||
|
A_neu = (1/tau) * M_neu + K_neu;
|
||||||
|
|
||||||
|
% Initial condition
|
||||||
|
u = u0;
|
||||||
|
|
||||||
|
for k = 1:Nt
|
||||||
|
b = (1/tau) * M_neu * u + F_neu;
|
||||||
|
u = A_neu \ b;
|
||||||
|
|
||||||
|
if mod(k,20) == 0
|
||||||
|
% figure(10)
|
||||||
|
% pdeplot(model, 'XYData', u, 'Mesh','on');
|
||||||
|
% axis equal
|
||||||
|
% title(['Insulated mug, t = ', num2str(k*tau), ' s']);
|
||||||
|
% colorbar
|
||||||
|
% drawnow
|
||||||
|
end
|
||||||
|
end
|
||||||
Loading…
Add table
Add a link
Reference in a new issue