Poisson Equation with Mixed Boundary Condition - MatlaB

[pinks]

Hi guys,
I'm solving a Poisson Equation with Mixed Boundary condition. But I have trouBle with that mixed BC in MATLAB. Anyone can help to fix? Thanks a lot!

dT^2/dx^2+dT/dy^2=-Q(x,y)/k

Rectangular domain (HxL), BC: Top: T(x,H)=Th, Left: dT/dx=0, Bottom: dT/dy=q, Right: dT/dx+B(T-Tinf)=0

Q(x,y) = 200*sin((2*pi/L)*x) + 175*sin((2*pi/H)*y)

nx ny L H q B k Th Tinf
20 30 2 3 -25 1 1 100 50

Coding: I tried 3 ways to define the mixed BC But it doesn't work.


Nx=20; %input('value of Nx:')
Ny=30; %input('value of Ny:')
L=2; %input('value of L:') %input initial value of Nx,Ny,L,H,Q,q,B
H=3; %input('value of H:')
q=-25; %input('value of q:')
B=1; %input('value of B:')
k=1; %input('value of k:')
Th=100; %input('value of Th:')
Tinf=50; %input('value of Tinf:')

err=1; tol=0.0001; iter=0; itmax=5000; %set initial error, tolerance, iteration limitation.
T=zeros(Nx+1,Ny+1);
T(:,Ny+1)=Th; % BOUNDARY CONDITION FOR TOP EDGE
dx=L/Nx; dy=H/Ny;
beta=(dx^2)/(dy^2);
x=0:dx:L;
y=0:dy;
Q=zeros(Nx+1,Ny+1);
for k=1:Nx+1
for l=1:Ny+1
Q(k,l)=200*sin((2*pi/L)*x(k)) + 175*sin((2*pi/H)*y(l));
end
end

while (err > tol)
iter=iter+1;
Told=T;
for i=1:Nx,j=1:Ny;
if i==Nx+1 ;
%---Tried 3 ways below but it doesn't work-----
%T(i,j) = (dx^2*Q(i,j)/k+(T(i-1,j)-2*dx*B*(T(i,j)-Tinf))+T(im1,j)+beta*T(i,jp1)+beta*T(i,jm1))/(2+2*beta);
%T(i,j) = (dx^2*Q(i,j)/k+2*T(i-1,j)+2*dx*B*Tinf+beta*T(i,jp1)+beta*T(i,jm1))/(2+2*beta+dx*B);
%T(i+1,j)=T(i-1,j)-2*dx*B*(T(i,j)-Tinf);
else ip1=i+1;
end
if i==1;
im1=2; % dT/dx=0 at left edge
else im1=i-1;
end
for j=1:Ny;
jp1=j+1;
if j==1;
jm1=2;
T(i,j) = (dx^2*Q(i,j)/k+ T(ip1,j)+T(im1,j)+beta*T(i,jp1)+beta*(T(i,jm1)-2*q*dy))/(2+2*beta);
else jm1=j-1;
end
T(i,j) = (dx^2*Q(i,j)/k+ T(ip1,j)+T(im1,j)+beta*T(i,jp1)+beta*T(i,jm1))/(2+2*beta);
%value of T we get from the equation
end
end
err=max(max(T-Told))/max(max(T));
if iter > itmax, break, end
end
Unu_Iterative=T';
[xp,yp]=meshgrid(x,y);
subplot(2,1,1);
surf(xp,yp,T')
subplot(2,1,2);
contour(x,y,T',20)

 

[NateD]

The code you have written is incomplete and needs to be fixed. The boundary condition for the right boundary should be implemented in your code. Also, the boundary condition for the bottom boundary is not correctly implemented. The boundary condition should be implemented as follows: T(i,j) = (dx^2*Q(i,j)/k + T(ip1,j)+T(im1,j)+beta*T(i,jp1) + beta*(T(i,jm1) - 2*q*dy))/(2+2*beta); T(Nx+1,j) = (dx^2*Q(Nx+1,j)/k + 2*T(Nx,j) + 2*dx*B*Tinf + beta*T(i,jp1) + beta*T(i,jm1))/(2+2*beta + dx*B);