Solving 2D Heat Equation on Non-Uniform Domain With Hole and Neumann Boundary Conditions
5 views (last 30 days)
Show older comments
Hi, I am having issues setting my Neumann Boundary Conditions for my problem. The domain is defined as
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1158263/image.png)
- 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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1158268/image.png)
With Boundary and Initial Conditions (Initial Condition take m = 0 inside the domain)
![On AB: ∂m/∂y=0
On BC: m=0.3
On CD: ∂m/∂n=0 where n is the normal direction
On DA: ∂m/∂x=0.
On HE: m=0.7
On GF: m=0.7
On EF: ∂m/∂y=0
On HG: ∂m/∂y=0
The initial conditions are to take m=0 inside the solution domain.](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1158273/image.png)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1158278/image.png)
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 :)
0 Comments
Answers (0)
See Also
Categories
Find more on Boundary Conditions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!