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

조회 수: 27(최근 30일)
Adam Harris
Adam Harris 2021년 11월 23일
댓글: VBBV 2021년 11월 25일
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
Adam Harris 2021년 11월 23일
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개)

VBBV
VBBV 2021년 11월 23일
편집: VBBV 2021년 11월 23일
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
VBBV
VBBV 2021년 11월 25일
Plus, some missing coefficient terms like (1-w) etc in Neumann BC

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by