When using contouf(X,Y,Z,_) with vectors X and Y and matrix Z of size m-by-n, the length of X must be n and the length of Y must be m. (Reference) In other words, the Y values correspond to the rows of Z, and the X values correspond to the columns of Z. Therefore, to fix your code, you can just transpose the matrix T when creating the contour plots.
% Parameters (units at m, C, W..)
L = 0.3; % Length of the cross section of square
H = 0.15; % Height of the cross section of square
T_flueinterior = 350; % Interior flue temperature
h_flueinterior = 100; % Interior flue convection coefficient
T_airexterior = 25; % Exterior air temperature
h_airexterior = 5; % Exterior air convection coefficient
k = 0.85; % Thermal conductivity of refractory brick
alpha = 5.5e-7; % Thermal diffusivity of refractory brick
Ti = 25; % Initial temperature of the wall
dx = 0.05; % Grid spacing at x direction
dy = 0.05; % Grid spacing at y direction
dt = 3600; % Time step in one hour in seconds
Nt = [1, 3, 10, 30]; % Time steps to output
% Converting Nt to seconds
Nt = Nt * 3600;
% Setting up the number of grid points in x and y direction
Nx = round(L/dx) + 1;
Ny = round(H/dy) + 1;
% Initializing temperature field
T = Ti * ones(Nx,Ny);
% Defining a coefficient to simplify the discretized FD equation
a = alpha * dt / dx^2;
% Time loop
for n = 1:max(Nt)/dt
T_new = T;
for j = 2:Ny-1
for i = 2:Nx-1
T_new(i,j) = T(i,j) + a * (T(i,j+1) + T(i,j-1) + T(i+1,j) + T(i-1,j) - 4 * T(i,j));
end
end
% Applying boundary conditions
% Left side flue
T_new(:,1) = T_flueinterior + (T(:,2) - T_flueinterior) / (1 + dx * h_flueinterior / k);
% Right side air
T_new(:,end) = T_airexterior + (T(:,end-1) - T_airexterior) / (1 + dx * h_airexterior / k);
% Top and bottom boundaries (adiabatic and insulated)
T_new(1,:) = T_new(2,:);
T_new(end, :) = T_new(end-1, :);
% Re-update temperature field
T = T_new;
% Save results at the specified times
if any(n*dt == Nt)
figure;
contourf(linspace(0,L,Nx), linspace(0,H,Ny), T.', 20, 'LineColor', 'none');
colorbar;
title(['Temperature Distribution at t = ', num2str(n*dt/3600), 'hours']);
xlabel('x (in m)');
ylabel('y (in m)');
end
end
% The final temperature distribution
figure;
contourf(linspace(0,L,Nx), linspace(0,H,Ny), T.', 20, 'LineColor', 'none');
title('Final Temperature Distribution');
xlabel('x (in m)');
ylabel('y (in m)');
