Plotting the 3D Heat Equation in 2D Slices

Hello, I am trying to plot the temperature distribution of a cake using the heat equation in 3D. The boundary conditions are below.

I am trying to plot the temperature distribution of a cake shaped as a rectangular prism on 2D slices. These slices are placed on the x-axis consecutively. (see below)

How can I visualize these slices? Is it also possible to make an animation of the temperature change over time for one of the slices? Thanks.

Answer 1

[Image Analyst]

"How can I visualize these slices? " <== There are several functions that can display a 2-D slice/matrix: imshow, imagesc, or image.

"Is it also possible to make an animation?" <==Once you've displayed them, you can get the currently displayed frame and save it as a frame in a movie. I'm attaching some demos. Look especially at movie_made_from_surf.m.

[uk]

Thanks for your response. Prior to your response, I have coded the heat equation to model the temperature distribution. Although it works, it is quite slow and when I try to increase the time step, I get stability issues and the output visuals seem weird. Do you know a way to fix them?

clear; clf;

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

alpha = rho / (cp * k);

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

Nx = 21; Nt = 35*6000;
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

[Torsten]

Your time step is too large. Look up "Courant condition" on how to set the maximum allowed time step size dt for Explicit Euler for the heat equation.

If I remember correctly, dt <= 0.25/(alpha*(1/dx^2+1/dy^2+1/dz^2)) where alpha is the Thermal Diffusivity.

By the way: Your alpha is incorrectly defined - it must read k/(rho*cp) instead of rho/(cp*k) (although for all values equal to 1 this doesn't matter).

rho = 1;
cp = 1;
k = 1;
alpha = k /( rho * cp);
Lx = 10;
Ly = 10;
Lz = 10;
Nx = 21; Nt = 35*6000;
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

dt = 0.0500

dtmax = 0.25/(alpha*(1/dx^2+1/dy^2+1/dz^2))

dtmax = 0.0208

[uk]

Changed the code as per your suggestions, thank you. The stability issues are gone. However, the issue still partly persists. It takes nearly 5 seconds for the simulation to evaluate 1 second. After a certain point, I can't really tell the difference either. Could the issue be about my implementation of the solution of the heat equation? Or is it related to my boundary conditions? Because I have realized at the beginning of the simulation the top color seems different than the other boundaries. For context, I am trying to run the simulation for about 35 simulated minutes. Thanks again.

Answer 2

[Torsten]

Your implementation of Explicit Euler looks correct.

If you want larger time steps and faster performance, use an ODE solver like ode15s.

[uk]

Isn't the heat equation a PDE though?

[Torsten]

Yes, but also in your code, you solve

dT(i,j,k)/dt = alpha* (((T(i+1,j,k) - 2*T(i,j,k) + T(i-1,j,k))/dx/dx)+...
                       ((T(i,j+1,k) - 2*T(i,j,k) + T(i,j-1,k))/dy/dy)+...
                       ((T(i,j,k+1) - 2*T(i,j,k) + T(i,j,k-1))/dz/dz))

thus (Nx-2)*(Ny-2)*(Nz-2) ODEs for the temperatures T(i,j,k) in the inner grid points, don't you ?

[uk]

I see. However, rather than solving the heat equation in ODE form, maybe I can use the PDE solver toolbox? It just seems easier to implement as I do not have much experience with MATLAB. Would there still be performance issues with this approach?

[Torsten]

I have no experience with the PDE toolbox, but I would assume that it's comparable in speed to an implementation in ODE15S.

[uk]

Alright, I'll try out PDE first. Thanks for helping me. One last (maybe not) question is could I plot the temperature distribution using the analytical solution I derived? What I mean is, I have derived the equation

where u_o and u_c is known. Maybe plotting the temperature using this equation could be a lot more efficient? How would I visualize 2D slices in this case?

[Torsten]

Evaluating a threefold series together with suitable stopping criteria to check when convergence has been achieved - I wouldn't go this way. Of course you can by setting x = 5, varying y and z in between 0 and 10 and setting t to the output time t = 35*60.