Neumann (zerogradient and open) boundary condition
11 views (last 30 days)
Show older comments
Hello all,
In my work of 2D rectangulat channel, right boundary is symmetry and top boundary is open to atmosphere. So, I need to use Neumann boundary condition.
At right side, u(:,Ny) = u(:,Ny-1);
At top side, u(1,:) = u(2,:);
At top and right boundary, the value is chosen from (N-1)th node.
The below code works well separately but when the same boundary condition is implemented in actual problem (involving numerical equations), I'm getting the velocities (u, v) as zeros at right and top boundaries.
Kindly someone share your ideas about it and guide me if I'm wrong anywhere.
Thank you
clear; clc; close all;
W =5; L=0.25;
Nx=20; %row
Ny=10; %column
x = linspace(0,L,Ny);
y = linspace(0,W,Nx);
[X, Y] = meshgrid(x,y);
dx = L/Ny; dy = W/Nx;
rhoo = 0.8; por = 0.3; mu =0.8E-5; dia = 0.003;
u = zeros(Nx, Ny);
v = zeros(Nx, Ny);
P = 101325 * ones(Nx, Ny);
GVecx = zeros(Nx,Ny);
GVecy = zeros(Nx,Ny);
%GVec = zeros(Nx,Ny);
% Boundary conditions
u(:,1) = 0; % left - 1st column @u=0
v(:,1) = 0; % left - 1st column @v=0
u(2,9) = 2;
u(3,9) = 3;
%for i = 2 : Nx-1
u(:,Ny) = u(:,Ny-1); % right - symmetry (Neumann condition)
v(:,Ny) = v(:,Ny-1); % right - symmetry (Neumann condition)
%end
%for j = 2 : Ny
u(1,:) = u(2,:); % top - open to atmosphere (Neumann condition)
v(1,:) = v(2,:); % top - open to atmosphere (Neumann condition)
%end
u(Nx,:) = 0; % bottom
v(Nx,:) = 0; % bottom
u(Nx,1) = 0; % bottom left corner(for left inlet condition @ u =0)
v(Nx,1) = 0; % bottom left corner(for left inlet condition @ u =0)
u(1,1) = 0; % top left corner (for left inlet condition @ v =0)
v(1,1) = 0; % top left corner (for left inlet condition @ v =0)
3 Comments
Torsten
on 18 Sep 2022
Edited: Torsten
on 18 Sep 2022
Maybe in your actual problem code, you initialize the velocities once at the start to 0, but don't update them as
u(:,Ny) = u(:,Ny-1); % right - symmetry (Neumann condition)
v(:,Ny) = v(:,Ny-1); % right - symmetry (Neumann condition)
%end
%for j = 2 : Ny
u(1,:) = u(2,:); % top - open to atmosphere (Neumann condition)
v(1,:) = v(2,:); % top - open to atmosphere (Neumann condition)
during the computation.
If you did update them, u(:,Ny-1), v(:,Ny-1), u(2,:) and v(2,:) also had to be 0 to make u(:,Ny), v(:,Ny), u(1,:) and v(1,:) equal 0 (which is not the case, I guess ?)
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!