Solving 2D Heat Equation on Non-Uniform Domain With Hole and Neumann Boundary Conditions

Hi, I am having issues setting my Neumann Boundary Conditions for my problem. The domain is defined as

• The lengths of the sides of the plate are AB = 1, AD = BC= 2/3.
• Point A is taken as the origin, and the coordinates of point E are (0.5,0.25).
• The lengths EF = HG = GF = HE = 0.3.
• The edge DC is given by the parabola

And the governing temperature equation is as such:

With Boundary and Initial Conditions (Initial Condition take m = 0 inside the domain)

In order to solve this problem, I want to transfer the coordinate system (x,y) to (), where this new system will yield a rectangular domain. Once this is done, I will use the Jacobi Method and Gauss Seidel to solve the problem. My current issue is with detting up the Neumann Boundary conditions. Attached are the funcitons I have written that pertain to the boundary conditions

function [X,Y] = boundary(N)

X = zeros(N); %% This function sets up the problem geometry.

Y = zeros(N);

%% External Boundary %%

y = @(t) (2/3) + (4/3)*(t-t.^2);

arclength = 1;%asinh(4/3)*(3/8) + 5/6; % arclength of given parabolic fcn

sum = 1 + 2*(2/3) + arclength; % sum of given lengths

AB = (1:round(N/sum));

BC = AB(end)+(1:round(N*(2/3)/sum));

CD = BC(end)+(1:round(N*arclength/sum));

DA = CD(end)+(1:round(N*(2/3)/sum));

v = size((1:round(N/sum)));

s = 1/v(1,end); % Sizing Factor

for i = AB

X(i,N) = (i-AB(1))*s;

Y(i,N) = 0;

end

for i = BC

X(i,N) = 1;

Y(i,N) = (i-BC(1))*s;

end

for i = CD

X(i,N) = 1 - (i-CD(1))*s;

Y(i,N) = y(X(i,N));

end

for i = DA

X(i,N) = 0;

Y(i,N) = (2/3) - (i-DA(1))*s;

end

Y(i,N) = 0;

%% Internal Boundary %%

E = [0.5,0.25];

length = 0.3;

diff = linspace(0,length,N/4);

zero = zeros(1,N/4);

X(:,1) = [E(1) + diff , E(1) + length + zero , E(1) + length - diff , E(1) + zero].';

Y(:,1) = [E(2) + zero , E(2) + diff, E(2) + length + zero , E(2) + length - diff].';

%% Combining and Plotting%%

for i=(1:N)

X(i,:) = linspace(X(i,1),X(i,N),N);

Y(i,:) = linspace(Y(i,1),Y(i,N),N);

end

figure,

hold on;

for i=(1:N)

plot(X(i,:),Y(i,:),'Color','r');

end

for j=(1:N)

plot(X([(1:N),1],j),Y([(1:N),1],j),'Color','b');

end

hold off

title('Problem Geometry')

xlabel('X Axis');

ylabel('Y Axis')

end

function T = conditions(T,N,X,Y) % This function develops the boundary and initial conditions.

N = 60; % I decided to define the boundary conditions on (X,Y) then later I will transfer coordinate systems.

T = zeros(N);

Ty = ones(N); % I defined Tx and Ty to simply see how the code responds to the boundary conditions, and made plots

Tx = ones(N); %The plots show that along the ecternal edges, it does not understand my boundary conditions.

arclength = 1;

sum = 1 + 2*(2/3) + arclength;

AB = (1:round(N/sum));

BC = AB(end)+(1:round(N*(2/3)/sum));

CD = BC(end)+(1:round(N*arclength/sum));

DA = CD(end)+(1:round(N*(2/3)/sum));

for i = AB

T(i,N) = T(i,N-1);

Ty(i,N) = (T(i,N) - T(i,N-1))/2;

Tx(i,N) = (T(N,i) - T(N-1,i))/2;

end

for i = BC

T(i,N) = -0.3;

Ty(i,N) = (T(i,N) - T(i,N-1))/2;

Tx(i,N) = (T(N,i) - T(N-1,i))/2;

end

for i = CD

nx = 8/3*X(i,N) - 4/3;

ny = 1;

T(i,N) = T(N-1,i)*ny + T(i,N-1)*nx;

Ty(i,N) = (T(i,N) - T(i,N-1))/2;

Tx(i,N) = (T(N,i) - T(N-1,i))/2;

end

for i = DA

T(N,i) = T(N-1,i);

Ty(N,i) = (T(i,N) - T(i,N-1))/2;

Tx(N,i) = (T(N,i) - T(N-1,i))/2;

end

diff = N/4;

for j = 1:diff

T(j,1) = T(j,N-1);

Ty(j,1) = (T(i,1) - T(i,N-1))/2;

Tx(j,1) = (T(1,i) - T(N-1,i))/2;

end

for j = j:diff*2

T(j,1) = 0.7;

Ty(j,1) = (T(i,1) - T(i,N-1))/2;

Tx(j,1) = (T(1,i) - T(N-1,i))/2;

end

for j = j:diff*3

T(j,1) = T(j,N-1);

Ty(j,1) = (T(i,1) - T(i,N-1))/2;

Tx(j,1) = (T(1,i) - T(N-1,i))/2;

end

for j = j:diff*4

T(j,1) = 0.7;

Ty(j,1) = (T(i,1) - T(i,N-1))/2;

Tx(j,1) = (T(1,i) - T(N-1,i))/2;

end

T(1,:) = (T(2,:) + T(N-1,:)) /2;

T(N,:) = T(1,:);

figure (1)

pcolor(X,Y,T)

colorbar;

figure (2)

pcolor(X,Y,Tx)

colorbar;

figure (3)

pcolor(X,Y,Ty)

colorbar;

end

Any help with this issue would be greatly appreciated!! Thank you :)