# How do I properly implement Neumann boundary condition for this problem?

87 views (last 30 days)
Adam Harris on 23 Nov 2021
Edited: Hamza SAYYOU on 24 Aug 2022
Hello world,
I've written a program to simulate channel flow past a rectangle using Successive Over-Relaxation. Everything is good except for my implimentation of a Neumann boundary condition.
The problem statement is:  My code for the upper left quadrant of the channel is as follows:
clear; close all; clc;
% Discretize The Domain & Initialize
h = 0.1; % Step size
x = 0:h:6; % x
y = 0:h:4; % y
Nx = length(x); % Number of steps in x
Ny = length(y); % Number of steps in y
PSI(Ny,Nx) = 0; % Initialize solution
% Define Acceleration Parameter for SOR
a = cos(pi/Nx) + cos(pi/Ny); % Alpha
w = (8-4*sqrt(4-a^2))/(a^2); % Optimal omega
%%% Define Boundary Conditions
for j=1:Ny
PSI(j,1) = y(j)/4; % PSI(y,0) Left end of channel, Dirichlet
end
for i=1:Nx
PSI(1,i) = 0; % PSI(0,x) Bottom of channel, Dirichlet
end
for j=1:Ny
PSI(j,Nx) = PSI(j,Nx-1); % PSI(y,6) Right end of channel **NEUMANN**
end
for i=1:Nx
PSI(Ny,i) = 1; % PSI(4,x) Top of channel, Dirichlet
end
%%% Implimentation of the Successive Over Relaxation Method
n=1000; % The number of iterations to be used
for i=1:n;
for ix = 2:Nx-1; % all internal x points
for iy = 2:Ny-1; % all internal y points
PSI(iy,ix) = (1.0-w)*PSI(iy,ix) + w/4.0*(PSI(iy,ix-1) + PSI(iy,ix+1) + PSI(iy-1,ix) + PSI(iy+1,ix));
if ix >= 51 && iy <= 21 % Boundary conditions of rectangle in channel
PSI(iy,ix) = 0;
end
end
end
end
Plotting the solution for the upper left quadrant of the channel, the issue is obvious on the right end of the domain. I've also solved the problem across the entire top half of the channel. This method avoids the Neumann BC due to axis-symmetry, so all of the BCs are Dirichlet. That solution plotted is: So I'm certain that the issue is my implimentation of the Neumann BC. Can someone please show me the proper way to incorporate the Neumann BC for this problem?
Thank you!

Adam Harris on 23 Nov 2021
I needed to plug PSI(j,Nx+1) = PSI(j,Nx-1) into the SOR equation at the index corresponding to the right end of the channel. So the SOR loop now looks like:
%%% Implimentation of the Successive Over Relaxation Method
for i=1:n;
for ix = 2:Nx-1; % all internal x points
for iy = 2:Ny-1; % all internal y points
PSI(iy,ix) = (1.0-w)*PSI(iy,ix) + w/4.0*(PSI(iy,ix-1) + PSI(iy,ix+1) + PSI(iy-1,ix) + PSI(iy+1,ix));
if ix >= 51 && iy <= 21 % Boundary conditions of rectangle in channel
PSI(iy,ix) = 0;
end
for j=2:Ny-1 % Neumann boundary condition on right end of the channel
PSI(j,Nx) = (1-w)*PSI(j,Nx) + (w/4)*(2*PSI(j,Nx-1) + PSI(j-1,Nx) + PSI(j+1,Nx));
end
end
end
end
and the solution plotted is now: ##### 1 CommentShowHide None
Hamza SAYYOU on 24 Aug 2022
Edited: Hamza SAYYOU on 24 Aug 2022
How do you plot the obstacle inside the contour ?

R2020a

### Community Treasure Hunt

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

Start Hunting!