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

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!

Answer 1

[Adam Harris]

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:

[Hamza SAYYOU]

Hello @Adam Harris,

How do you plot the obstacle inside the contour ?

[Hamza SAYYOU]

@Adam Harris can you please answer my question !!

[Torsten]

"please" and "!!" contradict each other :-)

[Shlok Vaibhav Singh]

Can you explain why PSI(j,Nx+1) = PSI(j,Nx-1) is preferred over PSI(j,Nx) = PSI(j,Nx-1) ? The former puts derivative at Nx=0 while the latter puts field at Nx-0.5 as 0, is that the reason?