MATLAB: Help with Understanding and Fixing Contour Plot

2d heat transfercontour plottingfinite element analysisMATLAB and Simulink Student Suite

Would anyone be able to help me understand how to fix this contour plot? I'm trying to get it to show the heat distribution for the attached object from a book problem. My code is below as well as the current output I'm getting from the contour. I've verified all the temperatures are correct for the temperature nodes. Any help would be appreciated.
clc;clear all;close all;
Lx = 5*10^-3;
Ly = 3*10^-3;
Nx = 6;
Ny = 4;
dx = Lx/Nx;
dy = Ly/Ny;
A = zeros(Nx*Ny);
B = zeros(Nx*Ny,1);
k = 25;
h_0 = 1000;
Tinf_0 = 1700;
h_i = 200;
Tinf_i = 400;
% Node 1
i = 1;
A(i,i) = -(2+(h_0*dx)/k);
A(i,i+1) = 1;
A(i,i+Nx) = 1;
B(i) = -((h_0*dx)/k)*Tinf_0;
% Nodes 2 - 5
for m = 2:Nx-1
    % n = 1
    i = m;
    A(i,i) = -2*(2+(h_0*dx)/k);
    A(i,i-1) = 1;
    A(i,i+1) = 1;
    A(i,i+Nx) = 2;
    B(i) = -2*((h_0*dx)/k)*Tinf_0;
end
% Node 6
i = Nx;
A(i,i) = -(2+(h_0*dx)/k);
A(i,i-1) = 1;
A(i,i+Nx) = 1;
B(i) = -((h_0*dx)/k)*Tinf_0;
% Nodes 7 & 13
for n = 2:Ny-1
    % m = 1 
    i = (n-1)*Nx + 1;
    A(i,i) = -4;
    A(i,i+1) = 2;
    A(i,i+Nx) = 1;
    A(i,i-Nx) = 1;
    B(i) = 0;
end
% Node 12
i = 12;
A(i,i) = -4;
A(i,i-1) = 2;
A(i,i+Nx) = 1;
A(i,i-Nx) = 1;
B(i) = 0;
% Node 20
i = 20;
A(i,i) = -4;
A(i,i+1) = 1;
A(i,i-1) = 1;
A(i,i-Nx) = 2;
B(i) = 0;
% Nodes 8 to 11
for m = 2:Nx-1
    for n = 2:Ny-2
        i = (n-1)*Nx + m;
        A(i,i) = -4;
        A(i,i+1) = 1;
        A(i,i-1) = 1;
        A(i,i+Nx) = 1;
        A(i,i-Nx) = 1;
        B(i) = 0;
    end
end
% Node 14
i = 14;
A(i,i) = -4;
A(i,i+1) = 1;
A(i,i-1) = 1;
A(i,i+Nx) = 1;
A(i,i-Nx) = 1;
B(i) = 0;
% Node 15
i = 15;
A(i,i) = -2*(3+(h_i*dx)/k);
A(i,i+1) = 1;
A(i,i-1) = 2;
A(i,i+Nx) = 1;
A(i,i-Nx) = 2;
B(i) = -2*((h_i*dx)/k)*Tinf_i;
% Nodes 16 & 17
for m = 4:Nx-1
    % n = 3
    i = (Ny-2)*Nx + m;
    A(i,i) = -2*(2+(h_i*dx)/k);
    A(i,i+1) = 1;
    A(i,i-1) = 1;
    A(i,i-Nx) = 2;
    B(i) = -2*((h_i*dx)/k)*Tinf_i;
end
% Node 18
i = 18;
A(i,i) = -(2+(h_i*dx)/k);
A(i,i-1) = 1;
A(i,i-Nx) = 1;
B(i) = -((h_i*dx)/k)*Tinf_i;
% Node 21
i = 21;
A(i,i) = -(2+(h_i*dx)/k);
A(i,i-1) = 1;
A(i,i-Nx) = 1;
B(i) = -((h_i*dx)/k)*Tinf_i;
% Node 19
i = 19;
A(i,i) = -2;
A(i,i+1) = 1;
A(i,i-Nx) = 1;
B(i) = 0;
% Nodes 22 to 24
for m = 4:Nx
    % n = Ny
    i = (Ny-1)*Nx + m; 
    A(i,i) = -Tinf_i;
    B(i) = -Tinf_i;
end
% disp(A)
% disp(B)
T = A\B;
T2d = zeros(4,6);
for n=1:Ny
    for m=1:Nx
        i = (n-1)*Nx + m;
        T2d(n,m) = T(i);
        x(m) = m*dx;
        y(n) = n*dy;
    end
end
X = x;
Y = y;
figure(1)
contourf(x,y,T2d,10)
% figure(2)
% plot(x,T(1:6))

Best Answer

I am not certain what you’re doing, or the result you want.
This sets the elements of ‘T2d’ that are equal to 1, to NaN, since NaN values do not plot, and that appears to be outside the area you want to plot.
See if that change, and figure(2) (the surf plot, rotated to match your ‘T2d’ matrix) do what you want. These lines replace your code after the ‘X’ and ‘Y’ assignments:
T2d(T2d == 1) = NaN;
figure(1)
contourf(x,y,T2d,10)
figure(2)
surf(T2d)
grid on
view(0, -90)
shading interp
Related Question