3D Heat Equation in MATLAB

Question

I am trying to model in MATLAB the temperature distribution inside a rectangular prism with boundary and initial conditions and heat equation I was trying to visualize 2D slices in the 3D shape. However, the simulation takes too long and when I try to change the time step, I get stability issues and the output visuals look weird and inconsistent. Is there a way to fix this? The code I used is below:

clear; clf;

rho = 1;
cp = 1;
k = 1;

alpha = rho / (cp * k);

Lx = 10;
Ly = 10;
Lz = 10;

Nx = 21; Nt = 35*60;
Ny = 21;
Nz = 21;
dx = Lx / (Nx - 1);
dy = Ly / (Ny - 1);
dz = Lz / (Nz - 1);

c = 1;
C = 0.1;
dt = C * dx / c;

Tn = zeros(Nx, Ny, Nz);
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
z = linspace(0, Lz, Nz);
[X, Y, Z] = meshgrid(x, y, z);

Tn(:, :, :) = 80;
t = 0;

Tn(1, :, :) = 350; Tn(end, :, :) = 350;
Tn(:, 1, :) = 350; Tn(:, end, :) = 350;
Tn(:, :, 1) = 350; Tn(:, :, end) = 350;

figure;
title(sprintf('Time = %f seconds', t));

for n = 1:Nt
    Tc = Tn;

    t = t + dt;
    for i = 2:Nx - 1
        for j = 2:Ny - 1
            for k = 2:Nz - 1
                Tn(i,j,k)=Tc(i,j,k) +...
                    dt * alpha *...
                    (((Tc(i+1,j,k) - 2*Tc(i,j,k) + Tc(i-1,j,k))/dx/dx)+...
                    ((Tc(i,j+1,k) - 2*Tc(i,j,k) + Tc(i,j-1,k))/dy/dy)+...
                    ((Tc(i,j,k+1) - 2*Tc(i,j,k) + Tc(i,j,k-1))/dz/dz));
            end
        end
    end

    % Display slices in the XZ plane
    slice(X, Y, Z, Tn, [], Lx/2, []);

    xlabel('X (inches)');
    ylabel('Y (inches)');
    zlabel('Z (inches)');
    title(sprintf('Time = %f seconds', t));

    % Pause for a short duration to visualize the changes
    pause(0.1);
end

I was trying to get an image like this. However, the simulation takes too long and after a short time there does not seem to be any temperature change whatsoever. When I increase the time step, I get a weird result. Do I have to wait for the simulation to take its time? And are there inconsistencies or flaws in my implementation of the heat equation in MATLAB?

Try first optimizing your code. for a 3D matrix `Tc`, `Tn(i,j,k)=Tc(i+1,j,k)-2*Tc(i,j,k)` is the same as (without loops) `Tn(1:end-1,:,:)=Tc(2:end,:,:)-2*Tc(1:end-1,:,:)`, for example. You can essentially remove the i,j,k loops to speed up your code significantly. — Ander Biguri, Jul 03 '23 at 10:26

Ander Biguri · Answer 1 · 2023-07-03T15:07:55.853

Several things with your code:

1- vectorize forms of the operation will speed up the code significantly. (see below)

2- you pause(0.1) to visualize. Nt*0.1=210, i.e. your code spends 3.5 minutes paused, therefore is slow. use drawnow.

3- Your dt is too big. If you run it with c=1 you will see it indeed goes weird. This is due to numerical inaccuracies, your simulation breaks. Reduce dt (e.g. c=5) and you will see this does not happen.

Code:

rho = 1;
cp = 1;
k = 1;

alpha = rho / (cp * k);

Lx = 10;
Ly = 10;
Lz = 10;

Nx = 21; Nt = 35*60;
Ny = 21;
Nz = 21;

dx = Lx / (Nx - 1);
dy = Ly / (Ny - 1);
dz = Lz / (Nz - 1);

c = 5;
C = 0.1;
dt = C * dx / c;

Tn = zeros(Nx, Ny, Nz);
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
z = linspace(0, Lz, Nz);
[X, Y, Z] = meshgrid(x, y, z);

Tn(:, :, :) = 80;
t = 0;

Tn(1, :, :) = 350; Tn(end, :, :) = 350;
Tn(:, 1, :) = 350; Tn(:, end, :) = 350;
Tn(:, :, 1) = 350; Tn(:, :, end) = 350;

figure;
title(sprintf('Time = %f seconds', t));

for n = 1:Nt
    Tc = Tn;

    t = t + dt;

    Tn=Tc;
    Tn(2:end-1,:,:)=Tn(2:end-1,:,:)+dt*alpha/dx^2*(Tc(3:end,:,:)-2*Tc(2:end-1,:,:)+Tc(1:end-2,:,:));
    Tn(:,2:end-1,:)=Tn(:,2:end-1,:)+dt*alpha/dy^2*(Tc(:,3:end,:)-2*Tc(:,2:end-1,:)+Tc(:,1:end-2,:));
    Tn(:,:,2:end-1)=Tn(:,:,2:end-1)+dt*alpha/dz^2*(Tc(:,:,3:end)-2*Tc(:,:,2:end-1)+Tc(:,:,1:end-2));

%     for i = 2:Nx - 1
%         for j = 2:Ny - 1
%             for k = 2:Nz - 1
%                 Tn(i,j,k)=Tc(i,j,k) +...
%                     dt * alpha *...
%                     (((Tc(i+1,j,k) - 2*Tc(i,j,k) + Tc(i-1,j,k))/dx/dx)+...
%                     ((Tc(i,j+1,k) - 2*Tc(i,j,k) + Tc(i,j-1,k))/dy/dy)+...
%                     ((Tc(i,j,k+1) - 2*Tc(i,j,k) + Tc(i,j,k-1))/dz/dz));
%             end
%         end
%     end

    % Display slices in the XZ plane
    h=slice(X, Y, Z, Tn, [], Lx/2, []);
    set(h,'edgecolor','none')    
    xlabel('X (inches)');
    ylabel('Y (inches)');
    zlabel('Z (inches)');
    title(sprintf('Time = %f seconds', t));

    % Pause for a short duration to visualize the changes
    drawnow
end

3D Heat Equation in MATLAB

1 Answers1