Correct location boundary conditions in heat transfer task

Not sure if this is the right place for this, but I'll try it anyway.

I have a code that enables me to perform a 2D heat transfer task on a concentric square plate (a square plate with a square hole in it). I use a Dirichlet boundary condition on the inner square and I face two issues (which might related).

The full code is shown below. Please run this to understand what I see.

close all
clear all

% Input
a = 2;
b = 1;
initial_temperature = 0;
inner_temperature = 50;
number_of_grids = 100;
amplitude_source = 0;
k = 1;

wide = a; % plate width
high = a; % plate height

resX = number_of_grids; % number of grid points in X
resY = number_of_grids; % number of grid points in Y

[X,Y] = meshgrid(0:wide/resX:wide , 0:high/resY:high); % Creates a rectangular grid

% Remove an inner square of bxb
mask = abs(X - (a/2)) <= (b/2) & abs(Y - (a/2)) <= (b/2);

X(mask) = NaN;
Y(mask) = NaN;

% Create a vector of the node positions
Xg = linspace(0,wide, resX+1);
Yg = linspace(0,high, resY+1);

% 2D-grid factors
dx = wide/resX; % grid size
dy = high/resY; % grid size
xweight = 1/dx^2; % weight factors
yweight = 1/dy^2; % weight factors
xyweight = 2/dx^2+2/dy^2;

% % Source term
Source = amplitude_source*sin(2*pi*X/wide).*sin(pi*Y/high); % This would be a sinusaoidally distributed source

% % Initial condition (all equal)
Phi(1:resY+1,1:resX+1) = initial_temperature;

iteration_multiplier = 100;
iterations = max(resX,resY);
relaxstep = linspace(1,iteration_multiplier*iterations,iteration_multiplier*iterations);

% Jacobi iteration
for i = 1:iteration_multiplier
    for j = 1:iterations
        Phitemp = Phi;
        % Interior points
        Phi(2:resY, 2:resX) = (yweight * (Phi(1:resY-1, 2:resX) + Phi(3:resY+1, 2:resX)) + ...
                               xweight * (Phi(2:resY, 1:resX-1) + Phi(2:resY, 3:resX+1)) + ...
                               Source(2:resY, 2:resX) / k) / xyweight;
        
        % Dirichlet on the inner square
        x1 = round((resX/2) - (((b/a)*resX)/2));
        x2 = round((resX/2) + (((b/a)*resX)/2));
        y1 = round((resY/2) - (((b/a)*resY)/2));
        y2 = round((resY/2) + (((b/a)*resY)/2));
        Phi(y1+1:y2+1, x1+1:x2+1) = inner_temperature;
    end
end

% Plot
figure
contourf(X,Y,Phi,20)
caxis([0 inner_temperature])
colorbar('horiz')
axis equal

The first issue is the placement of the boundary condition itself, which I have done inside the iteration. Currently I use

Phi(y1+1:y2+1, x1+1:x2+1) = inner_temperature;

but I am not sure if this is correct since I don't see the yellow (50) line when looking at the contourplot. When I change this to

Phi(y1:y2+2, x1:x2+2) = inner_temperature;

then I get to see the yellow line, but changing my a/b ratio at the input might cause this yellow line to be bigger or gone again. My questions are: Am I missing something here? Which way is correct and why?

The second issue has to do with the outer temperature. Currently it is zero, which would make sense since the initial is zero, but this is even when I change the a/b ratio to 9/8 (so a very slim shape). I believe this is quite a strange behaviour since it would imply that it doesn't matter if a plate is 1 meters wide or 1 cm wide: when I heat up one side it will always be zero on the other side of the plate. This is simply wrong, but I cannot find the thing in my code that causes this. I thought I was doing it correctly in the iteration with

Phi(2:resY, 2:resX) = (yweight * (Phi(1:resY-1, 2:resX) + Phi(3:resY+1, 2:resX)) + ...
                               xweight * (Phi(2:resY, 1:resX-1) + Phi(2:resY, 3:resX+1)) + ...
                               Source(2:resY, 2:resX) / k) / xyweight;

but it seems I am wrong or the result is correct but I don't see why.

Would be great if someone would know what I did wrong here.

Answer 1

[Karan Singh]

Hi Ganesh,

From what I understand you have provided a MATLAB code that simulates 2D heat transfer on a plate with a square hole in the center. In this simulation, you've encountered two main issues:

• Boundary Condition Placement: Your initial approach to set the inner temperature "(Phi(y1+1:y2+1, x1+1:x2+1) = inner_temperature;)" is correct. This applies the Dirichlet boundary condition to the nodes inside the inner square, excluding the edges. Changing it to "Phi(y1:y2+2, x1:x2+2) = inner_temperature; " incorrectly applies the condition to the edges. While for the yellow line disappearing I am not sure but won't at larger "a/b" ratios is that the heat transfer across the wider gap between the squares becomes more significant, leading to a smoother temperature distribution and a less pronounced boundary effect.
• Outer Temperature Behavior: Your current implementation correctly updates the temperature of interior nodes based on their neighboring nodes and the source term. However, it does not explicitly set the temperature of the outer boundary nodes. Consequently, these nodes retain their initial temperature (zero) throughout the iterations. To fix this, you need to define the desired temperature for the outer boundary nodes. This can be achieved by setting the elements of Phi corresponding to the edge nodes to the desired value (e.g., zero for a cold exterior or a specific temperature value) before starting the loop iterations. You can set them as

% Set outer temperature 
outer_temperature = 0; % Customize this value 
 
% Apply outer temperature on boundary nodes 
Phi(1,:) = outer_temperature; 
Phi(end,:) = outer_temperature; 
Phi(:,1) = outer_temperature; 
Phi(:,end) = outer_temperature;  

Hope this helps!

Karan Singh Khati