3D Heat Equation in PDE Solver

Hello. I have a code that solves my heat equation. However, when I interpolate temperatures at a specific point, I get an error:

Error using pde.ThermalResults.interpolateSolutionInternal

Incorrect number of input arguments.

Error in pde.TransientThermalResults/interpolateTemperature (line 208)

Tintrp = pde.ThermalResults.interpolateSolutionInternal(obj,varargin{:});

Error in untitled8 (line 42)

T_center = interpolateTemperature(results, 5,5,0.5);

I have the following code:

% Create a PDE model
model = createpde("thermal", "transient");

% Import the geometry from the STL file
g = importGeometry(model, "Plate10x10x1.stl");

% Plot the geometry with face labels and transparency
pdegplot(g, 'FaceLabels', 'on', 'FaceAlpha', 0.5);

% Customize the plot as needed
xlabel('x');
ylabel('y');
zlabel('z');
title('Cake Geometry');

% Define the boundary conditions
T_oven = 350;  % Oven temperature (in Fahrenheit)
T_initial = 80;  % Initial temperature of the cake (in Fahrenheit)

% Convert temperatures to Celsius for PDE Toolbox
T_oven = (T_oven - 32) * 5/9;
T_initial = (T_initial - 32) * 5/9;

% Define the thermal properties
thermalProperties(model, 'ThermalConductivity', 0.2, 'MassDensity', 1000, 'SpecificHeat', 4.2);

% Define the boundary conditions
thermalBC(model, 'Face', [1 2 3 4 5 6], 'Temperature', T_oven);

% Set the initial temperature distribution
thermalIC(model, T_initial);

% Generate the mesh
generateMesh(model);

tlist = 0:0.5:5;

% Solve the PDE
results = solve(model, tlist);

% Extract the temperature at the center of the cake
T_center = interpolateTemperature(results, 5,5,0.5);

% Convert the temperature back to Fahrenheit for plotting
T_center = T_center * 9/5 + 32;

% Plot the temperature of the center of the cake against time
t = results.Mesh.Nodes(4, :);
plot(t, T_center);
xlabel('Time');
ylabel('Temperature (°F)');
title('Temperature of the Center of the Cake');

% Convert the temperature to Celsius for display
T_center_celsius = (T_center - 32) * 5/9;
fprintf('Final temperature of the center of the cake: %.2f°C\n', T_center_celsius(end));

Answer 1

[Torsten]

T_center = interpolateTemperature(results, 5,5,0.5);

You didn't specify the times at which you want the temperature to be supplied. It's the fifth argument to "interpolateTemperature" that is missing.

Read the documentation for "interpolateTemperature" carefully:

https://uk.mathworks.com/help/pde/ug/pde.steadystatethermalresults.interpolatetemperature.html

[uk]

Thanks, I understand now.

[uk]

Hey Torsten, thanks for your generous helps so far, I appreciate it. Do you think the central temperature is supposed to be decreasing over time in this model? I was expecting it to increase over time.

[Torsten]

Do you think the central temperature is supposed to be decreasing over time in this model? I was expecting it to increase over time.

T_center should be >= T_initial for all x,y,z and t. But if T_center is slightly below T_initial, this can also happen due to numerical errors.

Thus I suggest you set rho and cp to 1 and lambda to a high value and see how temperature develops in the center point of your plate. It should approach the wall temperature quite quickly.

[uk]

What do you mean by lambda?

[Torsten]

lambda in the usual notation is Thermal Conductivity. You named it k.

[uk]

Fixed the code and it seems that the central temperature increases really fast. I don't think it should be this quick though. I guess it's about the dimensions of the "Plate10x10x1.stl" as I am comparing values based on different scales?

[Torsten]

I don't know. You might want to compare with your Euler-code. It was slow, but correct.

[uk]

Guess I'll stick with Euler method. I have one question though. The code I have written to show how the colors match temperatures in a legacy does not stay still. The upper limit does not seem to change but the value of the lower limit changes constantly throughout the simulation. How do I fix that?

% Constants
Lx = 9;  Ly = 13;  Lz = 2;               % Lengths of the cake in x, y, and z directions (in inches)
Nx = 91;  Ny = 131;  Nz = 41;          % Number of grid points in x, y, and z directions
T_initial = 80;                          % Initial temperature (in °F)
T_oven = 350;                            % Oven temperature (in °F)

% Create the grid
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
z = linspace(0, Lz, Nz);
[X, Y, Z] = meshgrid(x, y, z);

% Initialize the temperature array
T = T_initial * ones(Nx, Ny, Nz);
T(:, 1, :) = T_oven;  % Bottom boundary (oven temperature)
T(:, end, :) = T_oven;  % Top boundary (oven temperature)
T(:, :, 1) = T_oven;  % Front boundary (oven temperature)
T(:, :, end) = T_oven;  % Back boundary (oven temperature)
T(1, :, :) = T_oven;  % Left boundary (oven temperature)
T(end, :, :) = T_oven;  % Right boundary (oven temperature)

% Calculate the spatial steps
dx = Lx / (Nx - 1);
dy = Ly / (Ny - 1);
dz = Lz / (Nz - 1);

% Calculate the maximum stable time step based on CFL condition
dt_max = 1 / (2 * (1/dx^2 + 1/dy^2 + 1/dz^2));

% Set the time parameters
t_final = 45*60;         % Total simulation time (in seconds)
dt = min(dt_max, 1);       % Time step (in seconds), taking the smaller value between dt_max and 1
nt = floor(t_final / dt);

% Initialize the figure
figure;
xlabel('Y (inches)');
ylabel('Z (inches)');
title('Temperature Distribution (YZ Plane)');
axis equal;

% Set the color map for temperature
cmap = jet(4);  % Choose a colormap (adjust the number of colors as desired)
colormap(cmap);
clim([T_initial, T_oven]); % Set the color axis limits

% Create a meshgrid for the YZ plane
[y_mesh, z_mesh] = meshgrid(y, z);

% Plotting the initial temperature distribution
T_slice = squeeze(T(round(Nx/2), :, :));  % Extract the center slice (YZ plane)
h = imagesc(y, z, T_slice);
colorbar;
title('Temperature Distribution (YZ Plane)');

% Update the temperature distribution over time
for t = 1:nt
    T_new = T;
    for i = 2:Nx-1
        for j = 2:Ny-1
            for k = 2:Nz-1
                T_new(i, j, k) = T(i, j, k) + (dt / (dx^2)) * (T(i-1, j, k) + T(i+1, j, k) - 2*T(i, j, k)) + ...
                    (dt / (dy^2)) * (T(i, j-1, k) + T(i, j+1, k) - 2*T(i, j, k)) + ...
                    (dt / (dz^2)) * (T(i, j, k-1) + T(i, j, k+1) - 2*T(i, j, k));
            end
        end
    end
    T = T_new;

    % Extract the center slice (YZ plane)
    T_slice = squeeze(T(round(Nx/2), :, :));

    % Update the color map
    set(h, 'CData', T_slice);

    % Update the title
    title(sprintf('Temperature Distribution (YZ Plane) - Time: %.1f s', t));

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

[Torsten]

Use clim:

https://uk.mathworks.com/help/matlab/creating_plots/change-mapping-of-data-values-into-the-colormap.html