Solving 2D Heat Equation on Non-Uniform Domain With Hole and Neumann Boundary Conditions
5 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 :)
0 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Boundary Conditions finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!