steady state 2D Poisson's equation - boundary conditions

17 visualizaciones (últimos 30 días)
Anitha Limann
Anitha Limann el 7 de Mayo de 2023
Editada: Anitha Limann el 8 de Mayo de 2023
Hello,
I am trying to solve steady state 2D Poisson's equation for a two spreading ridges problem. I have introduced boundary conditions tto. However though I should not have a flow across the top and bottom boundaries as you can see in the resultant image i do get flow across those boundaries. This gives me bulls eye contours on the spreading ridges but I should be getting horse shoe shapes instead.
This code solves for the 2D steady-state Poisson’s equation using finite difference method. The Poisson’s equation is given by:
∇²phi(x,y)=Gstar(x,y)
where phi(x,y) is the unknown function and Gstar(x,y) is a known function. The solution is obtained on a square domain with Dirichlet boundary conditions. The potential should not have this bull's eye shape. Potential about the ridges should be low and increase towards the boundaries. At the top and the bottom boundaries near ridges I know I should get a shape similar to horse shoe but as you can see something is not correct in my boundary conditions. For the top and bottom boundaries dphi/dx should be 0, which I thought I am doing but I cannot get this to work.
I am not sure if this is a math problem or Matlab problem anymore. Please help me with solving this problem.
%% Define grid parameters
x = 0:5:500; dx=x(2)-x(1); nx=length(x);
y = 0:5:500; dy=y(2)-y(1); ny=length(y);
% Define ridge locations
r1_x = 200;
r1_y = 250:500;
r2_x = 400;
r2_y = 0:250;
% Create meshgrid for plotting
[xx,yy] = meshgrid(x,y);
upper_half = yy >= 250;
lower_half = yy < 250;
% Apply age values to upper and lower halves of grid
t_upper = abs(xx - r1_x)/4;
t_lower = abs(xx - r2_x)/4;
t = zeros(size(xx));
t(upper_half) = t_upper(upper_half);
t(lower_half) = t_lower(lower_half);
t=max(t,0.1);
% Compute derivative of g
G = 0.11./(2*sqrt(t));
% Ignore infinity VALUES
count=0;
Gsum=0;
for i=1:nx
for j=1:ny
if G(i,j)<=0.11
Gsum=Gsum+G(i,j);
count=count+1;
end
end
end
% Calculate the mean
Gmean=Gsum/count;
Gstar=G-Gmean;
%% Compute potentials using Poisson's equation (del^2*potential = Gstar)
phi = zeros(nx,ny); % Initialize solution matrix
% Top boundary condition
for j=2:nx-1
phi(1,j) =phi(1,j-1)+phi(1,j+1)+phi(2,j)+phi(2,j)-(dx^2*Gstar(1,j)/4);
end
% Bottom boundary condition
for j=2:nx-1
phi(ny,j) =phi(ny,j-1)+phi(ny,j+1)+phi(ny-1,j)+phi(ny-1,j)-(dx^2*Gstar(ny,j)/4);
end
% left boundary condition
for i=2:ny-1
phi(i,1) =phi(i,2)+phi(i,2)+phi(i-1,1)+phi(i+1,1)-(dy^2*Gstar(i,1)/4);
end
% right boundary condition
for i=2:ny-1
phi(i,end) =phi(i,end-1)+phi(i,end-1)+phi(i-1,end)+phi(i+1,end)-(dy^2*Gstar(i,end)/4);
end
% Top right corner
phi(1,1) = (dx^2*dy^2/(2*(dx^2 + dy^2))) * ((2*(phi(2,1))/dx^2)+ (2*(phi(1,2))/dy^2) - Gstar(1,1)) ;
% Top left corners
phi(1,nx) = (dx^2*dy^2/(2*(dx^2 + dy^2))) * ((2*(phi(2,nx))/dx^2)+ (2*(phi(1,nx-1))/dy^2) - Gstar(1,nx)) ;
% Bottom right corner
phi(ny,1) = (dx^2*dy^2/(2*(dx^2 + dy^2))) * ((2*(phi(ny-1,1))/dx^2)+ (2*(phi(ny,2))/dy^2) - Gstar(ny,1));
% Bottom left corner
phi(ny,nx) = (dx^2*dy^2/(2*(dx^2 + dy^2))) * ((2*(phi(ny-1,nx))/dx^2)+ (2*(phi(ny,nx-1))/dy^2) - Gstar(ny,nx));
%%
% Jacobi iteration method
maxiter = 10000; % maximum number of iterations
tol = 1e-6; % tolerance for convergence
for iter = 1:maxiter
phi_old = phi;
for i = 2:ny-1
for j = 2:nx-1
phi(i,j) = (dx^2*dy^2/(2*(dx^2 + dy^2))) * (((phi(i+1,j) + phi(i-1,j))/dx^2)+ ((phi(i,j+1) + phi(i,j-1))/dy^2) - Gstar(i,j)) ;
end
end
E = max(abs(phi(:)-phi_old(:)));
if E < tol % convergence check
break;
end
end
%% Convert potential into velocities (u and v components)
p = reshape(phi,ny,nx);
u = 0*p;
v = 0*p;
u(:,2:end-1) = -(p(:,3:end) - p(:,1:end-2))/(2*dx);
v(2:end-1,:) = (p(3:end,:) - p(1:end-2,:))/(2*dy);
scale=5;
figure; hold on;axis fill;
q=quiver(x,y,u,v,scale);
contour(x,y,p,LineWidth=1)
  3 comentarios
Anitha Limann
Anitha Limann el 7 de Mayo de 2023
Thank you for your kind answer. When i include boundary conditions in to jacobi interation method I get empty plots.
Torsten
Torsten el 7 de Mayo de 2023
Editada: Torsten el 8 de Mayo de 2023
When i include boundary conditions in to jacobi interation method I get empty plots.
Most probably your Jacobi iterations diverge and you get NaN values for p which then are not plotted. But getting no plot is not a reason not to code correctly.
And what about changing phi to phi_old for the right-hand side ?
And what about a problem description ?

Iniciar sesión para comentar.

Respuestas (0)

Categorías

Más información sobre Loops and Conditional Statements en Help Center y File Exchange.

Productos

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by