Correct location boundary conditions in heat transfer task

3 次查看（过去 30 天）
Ganesh Gebhard
Ganesh Gebhard 2023-11-12
回答： Karan Singh 2023-12-13
Ran in:
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.

请先登录，再回答此问题。

采纳的回答

Karan Singh
Karan Singh 2023-12-13
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by