Solving 2D Heat Equation on Non-Uniform Domain With Hole and Neumann Boundary Conditions
조회 수: 5 (최근 30일)
이전 댓글 표시
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
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Boundary Conditions에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!