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

28 views (last 30 days)
Adam Harris
Adam Harris on 23 Nov 2021
Commented: VBBV on 25 Nov 2021
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!

Accepted Answer

Adam Harris
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:

More Answers (1)

VBBV
VBBV on 23 Nov 2021
Edited: VBBV on 23 Nov 2021
for j=1:Ny-1
PSI(j,Nx) = (PSI(j+1,Nx-1)- PSI(j,Nx-1))/(x(Nx)-x(Nx-1)); % PSI(y,6) Right end of channel **NEUMANN**
end
Try this and check
  6 Comments

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!