Why my heat node goes below zero??

% Simulation Parameters
dt = 10; % Time step [s]
t_max = 2000; % Maximum simulation time [s]
Nx = 34; % Number of horizontal nodes (가로)
Ny = 40; % Number of vertical nodes (세로)
dx = 1; dy = 1; % Spatial resolution
T_infinite = 20; % Ambient air temperature
sigma = 5.67e-8; % Stefan-Boltzmann constant
epsilon = 0.3;
h_air = 15; % Convective heat transfer coefficient [W/m^2K]
h = ones(Ny, Nx); % Initialize convective heat transfer coefficient
% Initialize material properties
k = 0.025 * ones(Ny, Nx); % Default thermal conductivity
c = ones(Ny, Nx); % Specific heat capacity
rho = 1.2 * ones(Ny, Nx);
T = T_infinite * ones(Ny, Nx); % Initial temperature field
T(i, end)=400;
T_conv_bottom=20;
T_conv_left=20;
% ABS (x=6:8, y=1:13), (x=28:30, y=1:13), (x=9:27, y=1:4)
k(1:13, 6:8) = 0.23; T(1:13, 6:8) = 21; rho(1:13, 6:8) = 1040;
k(1:13, 28:30) = 0.23; T(1:13, 28:30) = 21; rho(1:13, 28:30) = 1040;
k(1:4, 9:27) = 0.23; T(1:4, 9:27) = 21; rho(1:4, 9:27) = 1040;
c(1:13, 6:8) = 1.5; % 비열 [J/gK]
c(1:13, 28:30) = 1.5;
c(1:4, 9:27) = 1.5;
% Handle (x=6:9, y=29:30), (x=6:7, y=20:28)
k(29:30, 6:9) = 0.3; T(29:30, 6:9) = 19; rho(29:30, 6:9) = 950;
k(20:28, 6:7) = 0.3; T(20:28, 6:7) = 19; rho(20:28, 6:7) = 950;
c(29:30, 6:9) = 2.0;
c(20:28, 6:7) = 2.0;
% Vacuum (x=11, y=10:31), etc.
h(10:31, 11) = 0.0003; 
h(10, 11:12) = 0.0003; 
h(6:10, 12) = 0.0003; 
h(6, 12:24) = 0.0003; 
h(6:10, 24) = 0.0003; 
h(10, 24:25) = 0.0003; 
T(10:31, 25) = 14; 
T(10:31, 11) = 14; 
T(10, 11:12) = 14; 
T(6:10, 12) = 14; 
T(6, 12:24) = 14; 
T(6:10, 24) = 14; 
T(10, 24:25) = 14; 
T(10:31, 25) = 14;
rho(10:31, 11) = 0;
rho(10, 11:12) = 0;
rho(6:10, 12) = 0;
rho(6, 12:24) = 0;
rho(6:10, 24) = 0;
rho(10, 24:25) = 0;
c(10:31, 11) = 0;
c(10, 11:12) = 0;
c(6:10, 12) = 0;
c(6, 12:24) = 0;
c(6:10, 24) = 0;
c(10, 24:25) = 0;
% Ice (x=15:21, y=27:30)
k(27:30, 15:21) = 2.2; c(27:30, 15:21) = 2.06; T(27:30, 15:21) = -15;
rho(27:30, 15:21) = 917;
c(27:30, 15:21) = 2.06;
% Water (x=13, y=12:30), etc.
k(12:30, 13) = 0.575; h(12:30, 13) = 200; c(12:30, 13) = 4.186; T(12:30, 13) = 4;
k(8:30, 14) = 0.575; h(8:30, 14) = 200; c(8:30, 14) = 4.186; T(8:30, 14) = 4;
k(8:26, 15:21) = 0.575; h(8:26, 15:21) = 200; c(8:26, 15:21) = 4.186; T(8:27, 15:21) = 4;
k(8:30, 22) = 0.575; h(8:30, 22) = 200; c(8:30, 22) = 4.186; T(8:30, 22) = 4;
k(12:30, 23) = 0.575; h(12:30, 23) = 200; c(12:30, 23) = 4.186; T(12:30, 23) = 4;
rho(12:30, 13) = 1000;
rho(8:30, 14) = 1000;
rho(8:36, 15:21) = 1000;
rho(8:30, 22) = 1000;
c(12:30, 23) = 4.186;
c(12:30, 13) = 4.186;
c(8:30, 14) = 4.186;
c(8:36, 15:21) =4.186;
c(8:30, 22) = 4.186;
c(12:30, 23) =4.186;
% Air inside tumbler (x=15:21, y=31:34), etc.
k(31:34, 15:21) = 0.025; h(31:34, 15:21) = 2;
k(31:32, 14) = 0.025; h(31:32, 14) = 2;
k(31, 13) = 0.025; h(31, 13) = 2;
k(31, 23) = 0.025; h(31, 23) = 2;
k(31:32, 22) = 0.025; h(31:32, 22) = 2;
T(31:32, 14) = 8;
T(31, 13) = 8;
T(31, 23) = 8;
T(31:32, 22) = 8; 
T(31:34,15:21) = 8;
rho(31:34, 15:21) = 1.2;
rho(31:32, 14) = 1.2;
rho(31, 13) = 1.2;
rho(31, 23) = 1.2;
rho(31:32, 22) = 1.2;
c(31:34, 15:21) = 1.005;
c(31:32, 14) = 1.005;
c(31, 13) = 1.005;
c(31, 23) = 1.005;
c(31:32, 22) = 1.005;
% Stainless steel outside (x=12:14, y=34), etc.
k(34, 12:14) = 15.2; T(34, 12:14) = 19;
k(34, 22:24) = 15.2; T(34, 22:24) = 19;
k(33, 11:12) = 15.2; T(33, 11:12) = 19;
k(33, 24:25) = 15.2; T(33, 24:25) = 19;
k(32, 25:26) = 15.2; T(32, 25:26) = 19;
k(32, 10:11) = 15.2; T(32, 10:11) = 19;
k(9:32, 26) = 15.2; T(9:32, 26) = 19;
k(31:32, 11) = 15.2; T(31:32, 10) = 19;
k(9:32, 10) = 15.2; T(9:32, 10) = 19;
k(5:9, 11) = 15.2; T(5:9, 11) = 19;
k(5:9, 25) = 15.2; T(5:9, 25) = 19;
k(5, 11:25) = 15.2; T(5, 11:25) = 19;
rho(34, 12:14) = 8000; 
rho(34, 22:24) = 8000; 
rho(33, 11:12) = 8000; 
rho(33, 24:25) = 8000; 
rho(32, 25:26) = 8000;
rho(32, 10:11) = 8000;
rho(9:32, 26) = 8000; 
rho(31:32, 11) = 8000;
rho(9:32, 10) = 8000; 
rho(5:9, 11) = 8000;
rho(5:9, 25) = 8000; 
rho(5, 11:25) = 8000;
c(34, 12:14) = 0.5; 
c(34, 22:24) = 0.5; 
c(33, 11:12) = 0.5; 
c(33, 24:25) = 0.5; 
c(32, 25:26) =0.5;
c(32, 10:11) = 0.5;
c(9:32, 26) = 0.5; 
c(31:32, 11) =0.5;
c(9:32, 10) = 0.5; 
c(5:9, 11) = 0.5;
c(5:9, 25) = 0.5; 
c(5, 11:25) = 0.5;
% Stainless steel inside (x=13:14, y=33), etc.
k(33, 13:14) = 16.2; T(33, 13:14) = 5;
k(33, 22:23) = 16.2; T(33, 22:23) = 5;
k(32, 23:24) = 16.2; T(32, 23:24) = 5;
k(32, 12:13) = 16.2; T(32, 12:13) = 5;
k(11:31, 12) = 16.2; T(11:31, 12) = 5;
k(11:31, 24) = 16.2; T(11:31, 24) = 5;
k(11, 23:24) = 16.2; T(11, 23:24) = 5;
k(11, 12:13) = 16.2; T(11, 12:13) = 5;
k(7:10, 13) = 16.2; T(7:10, 13) = 5;
k(7:10, 23) = 16.2; T(7:10, 23) = 5;
k(7, 13:23) = 16.2; T(7, 13:23) = 5;
rho(33, 13:14) = 8000; 
rho(33, 22:23) = 8000; 
rho(32, 23:24) = 8000; 
rho(32, 12:13) = 8000; 
rho(11:31, 12) = 8000; 
rho(11:31, 24) = 8000;
rho(11, 23:24) = 8000; 
rho(11, 12:13) = 8000;
rho(7:10, 13) = 8000; 
rho(7:10, 23) = 8000;
rho(7, 13:23) = 8000; 
c(33, 13:14) = 0.5; 
c(33, 22:23) = 0.5; 
c(32, 23:24) = 0.5; 
c(32, 12:13) = 0.5; 
c(11:31, 12) = 0.5; 
c(11:31, 24) = 0.5; 
c(11, 23:24) = 0.5; 
c(11, 12:13) = 0.5; 
c(7:10, 13) = 0.5; 
c(7:10, 23) = 0.5; 
c(7, 13:23) = 0.5; 
%lid
k(35:36,14:22)= 0.2;
T(35:36,14:22)= 15;
rho(35:36, 14:22) = 950;
c(35:36, 14:22) = 1.75;
% Simulation Loop
figure;
colormap(jet);
imagesc(flipud(T), [-20, 400]);
colorbar;
title('Initial Temperature Distribution (t = 0)');
xlabel('X-Axis (Nodes)');
ylabel('Y-Axis (Nodes)');
videoFile = 'simulation_result.mp4';
v = VideoWriter(videoFile, 'Motion JPEG AVI');
v.FrameRate = 10;
open(v);
for t = 1:t_max
    T_new = T; % Create a new array for updated temperatures
    
    % Update internal nodes
    for i = 2:Ny-1
        for j = 2:Nx-1
            T_cond = k(i, j) / (rho(i, j) * c(i, j)) * ...
                ((T(i+1, j) - 2*T(i, j) + T(i-1, j)) / dx^2 + ...
                (T(i, j+1) - 2*T(i, j) + T(i, j-1)) / dy^2);
            
            T_conv = h(i, j) * (T_infinite - T(i, j)) / (rho(i, j) * c(i, j));
            
            T_new(i, j) = T(i, j) + dt * (T_cond + T_conv);
        end
    end
    
    % Boundary conditions
    % Bottom boundary (convection)
    for j = 1:Nx
        T_conv_bottom = h_air / (rho(Ny, j) * c(Ny, j)) * (T_infinite - T(Ny, j));
        T_new(Ny, j) = T(Ny, j) + dt * T_conv_bottom;
    end
    
    % Left boundary (convection)
    for i = 1:Ny
        T_conv_left = h_air / (rho(i, 1) * c(i, 1)) * (T_infinite - T(i, 1));
        T_new(i, 1) = T(i, 1) + dt * T_conv_left;
    end
    
    % Top boundary (insulated)
    T_new(1, 2:Nx-1) = T(2, 2:Nx-1);
    % Right boundary (fixed 400°C and radiation from 38 to 37)
    T_new(:, Nx) = 400;
    for i = 2:Ny-1
        if i == 38
            q_rad = epsilon * sigma * ((T(i, Nx) + 273.15)^4 - (T(i, Nx-1) + 273.15)^4);
            T_new(i, Nx-1) = T_new(i, Nx-1) + dt * q_rad / (rho(i, Nx-1) * c(i, Nx-1));
        end
    end
    % Update temperature field
    T = T_new;
    % Visualization
    if mod(t, 10) == 0
        imagesc(flipud(T), [-20, 30]);
        colorbar;
        title(['Time = ', num2str(t * dt), ' seconds']);
        xlabel('X-Axis (Nodes)');
        ylabel('Y-Axis (Nodes)');
        drawnow;
        frame = getframe(gcf);
        writeVideo(v, frame);
    end
end
close(v);

[Torsten]

I suggest using the PDE Toolbox.

[Sandeep Mishra]

Hi @준 최

Could you please specify the value or range of i for the line T(i, end) = 400; ?

This will help in diagnosing the error. Thank you!

[준 최]

i changed the code and reuploaded whats the matter

[Sandeep Mishra]

Hi @준 최

Can you update the existing code or share the variable for easier diagnosis? Also, can you specify the issue you're encountering and your expected result?