# Open Boundary Condition

Dear all,
I have problem with open boundary condition implementation.
I have changed the reference MATLAB code in palabos for flow in a channel. but when I want to implement open BC (East boundary condition in code below, fIn(7),fIn(4), fIn(8)), the result becomes unstable.
Any help will be appreciated kindly.
My code is as follows:

``````
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flow_in_Channel.m: 2D channel flow, simulated by a LB method           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lattice Boltzmann sample, Matlab script                                %
% Implementation of 2d channel geometry and Zou/He boundary              %
% condition by Emad Tahmoures                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Purpose:                                                               %
%   This program simulates the flow in a channel by                      %
%   Lattice-Boltzmann method for a 2D Channel problem.                   %
%                                                                        %
% Lattice structure:                                                     %
%  c7  c3   c6                                                           %
%    \  |  /                                                             %
%  c4 -c1 - c2                                                           %
%    /  |  \                                                             %
%  c8  c5   c9                                                           %
%                                                                        %
% Record of revisions:                                                   %
%                                                                        %
%      Date           Programmer            Description of change        %
%     ======         ============          =======================       %
%  15 March 2014   Adriano Sciacovelli          Original code            %
%                                                                        %
% Last Revision: 15 March 2014                                           %
%                                                                        %
% Define variables:                                                      %
%   lx          -- Number of lattice nodes in x direction                %
%   ly          -- Number of lattice nodes in y direction                %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc; clear; tic;

% GENERAL FLOW CONSTANTS
lx = 100;
ly = 40;

uChannel  = 0.05; % horizontal channel velocity
vChannel  = 0;    % vertical channel velocity
Re    = 100;  % Reynolds number
nu    = uChannel *lx / Re;     % kinematic viscosity
omega = 1. / (3*nu+1./2.); % relaxation parameter
maxT  = 400; % total number of iterations
tPlot = 10;    % cycles for graphical output

% D2Q9 LATTICE CONSTANTS
t   = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cx  = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
cy  = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
opp = [ 1, 4, 5, 2,  3, 8, 9,  6,  7];
channel = [2:(ly-1)]; % fluid region

[y,x] = meshgrid(1:ly,1:lx);
obst = ones(lx,ly);
obst(1:lx,channel) = 0;
bbRegion = find(obst); % boundary of the lid

% INITIAL CONDITION: (rho=0, u=0) ==> fIn(i) = t(i)
fIn = reshape( t' * ones(1,lx*ly), 9, lx, ly);

% MAIN LOOP (TIME CYCLES)
for cycle = 1:maxT

% MACROSCOPIC VARIABLES
rho = sum(fIn);
ux = reshape ( (cx * reshape(fIn,9,lx*ly)), 1,lx,ly ) ./rho;
uy = reshape ( (cy * reshape(fIn,9,lx*ly)), 1,lx,ly ) ./rho;

% MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS

ux(:,1,channel) = uChannel; % channel inlet x-velocity
uy(:,1,channel) = vChannel; % channel inlet y-velocity
rho(:,1,channel) = 1 ./ (1-ux(:,1,channel)) .* ( ...
sum(fIn([1,3,5],1,channel)) + 2*sum(fIn([4,7,8],1,channel)) );

% MICROSCOPIC BOUNDARY CONDITIONS: West boundary condition (Zou/He BC)
fIn(2,1,channel) = fIn(4,1,channel) + 2/3*rho(:,1,channel).*ux(:,1,channel);
fIn(6,1,channel) = fIn(8,1,channel) + 1/2*(fIn(5,1,channel)-fIn(3,1,channel))+ ...
1/6*rho(:,1,channel).*ux(:,1,channel) + 1/2*rho(:,1,channel).*uy(:,1,channel);
fIn(9,1,channel) = fIn(7,1,channel) + 1/2*(fIn(3,1,channel)-fIn(5,1,channel))+ ...
1/6*rho(:,1,channel).*ux(:,1,channel) - 1/2*rho(:,1,channel).*uy(:,1,channel);

% % MICROSCOPIC BOUNDARY CONDITIONS: East boundary condition
fIn(7,lx,channel)=2*fIn(7,lx-1,channel)-fIn(7,lx-2,channel);
fIn(4,lx,channel)=2*fIn(4,lx-1,channel)-fIn(4,lx-2,channel);
fIn(8,lx,channel)=2*fIn(8,lx-1,channel)-fIn(8,lx-2,channel);

% COLLISION STEP
for i=1:9
cu = 3*(cx(i)*ux+cy(i)*uy);
fEq(i,:,:) = rho .* t(i) .* ...
( 1 + cu + 1/2*(cu.*cu) - 3/2*(ux.^2+uy.^2) );
fOut(i,:,:) = fIn(i,:,:) - omega .* (fIn(i,:,:)-fEq(i,:,:));
end

% MICROSCOPIC BOUNDARY CONDITIONS: NO-SLIP WALLS (bounce-back)
for i=1:9
fOut(i,bbRegion) = fIn(opp(i),bbRegion);
end

% STREAMING STEP
for i=1:9
fIn(i,:,: ) = circshift(fOut(i,:,: ), [0,cx(i),cy(i)]);
end

% VISUALIZATION
if (mod(cycle,tPlot)==0)
u = reshape(sqrt(ux.^2+uy.^2),lx,ly);
u(bbRegion) = nan;
imagesc(u(:,ly:-1:1)'./uChannel);
colorbar
axis equal off; drawnow
end

end
toc;

``````

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flow_in_Channel.m: 2D channel flow, simulated by a LB method %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lattice Boltzmann sample, Matlab script %
% Implementation of 2d channel geometry and Zou/He boundary %
% condition by Emad Tahmoures %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Purpose: %
% This program simulates the flow in a channel by %
% Lattice-Boltzmann method for a 2D Channel problem. %
% %
% Lattice structure: %
% c7 c3 c6 %
% \ | / %
% c4 -c1 - c2 %
% / | \ %
% c8 c5 c9 %
% %
% Record of revisions: %
% %
% Date Programmer Description of change %
% ====== ============ ======================= %
% 15 March 2014 Adriano Sciacovelli Original code %
% %
% Last Revision: 15 March 2014 %
% %
% Define variables: %
% lx – Number of lattice nodes in x direction %
% ly – Number of lattice nodes in y direction %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc; clear; tic;

% GENERAL FLOW CONSTANTS
lx = 800;
ly = 40;

uChannel = 0.05; % horizontal channel velocity
vChannel = 0; % vertical channel velocity
Re = 100; % Reynolds number
nu = uChannel lx / Re; % kinematic viscosity
omega = 1. / (3
nu+1./2.); % relaxation parameter
maxT = 40000; % total number of iterations
tPlot = 10; % cycles for graphical output

% D2Q9 LATTICE CONSTANTS
t = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cx = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
cy = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
opp = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];
channel = (2:(ly-1)); % fluid region

[y,x] = meshgrid(1:ly,1:lx);
obst = ones(lx,ly);
obst(1:lx,channel) = 0;
bbRegion = find(obst); % boundary of the lid

% INITIAL CONDITION: (rho=0, u=0) ==> fIn(i) = t(i)
fIn = reshape( t’ * ones(1,lx*ly), 9, lx, ly);
fOut=zeros(9,lx,ly);
fEq=zeros(9,lx,ly);
% MAIN LOOP (TIME CYCLES)
for cycle = 1:maxT

% MACROSCOPIC VARIABLES
rho = sum(fIn);
ux = reshape ( (cx * reshape(fIn,9,lxly)), 1,lx,ly ) ./rho;
uy = reshape ( (cy * reshape(fIn,9,lx
ly)), 1,lx,ly ) ./rho;

% MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS

ux(:,1,channel) = uChannel; % channel inlet x-velocity
uy(:,1,channel) = vChannel; % channel inlet y-velocity
rho(:,1,channel) = 1 ./ (1-ux(:,1,channel)) .* ( …
sum(fIn([1,3,5],1,channel)) + 2*sum(fIn([4,7,8],1,channel)) );

% COLLISION STEP
for i=1:9
cu = 3*(cx(i)ux+cy(i)uy);
fEq(i,:, = rho .
t(i) .

( 1 + cu + 1/2*(cu.cu) - 3/2(ux.^2+uy.^2) );
fOut(i,:, = fIn(i,:, - omega .* (fIn(i,:,:)-fEq(i,:,:));
end

% MICROSCOPIC BOUNDARY CONDITIONS: NO-SLIP WALLS (bounce-back)
for i=1:9
fOut(i,bbRegion) = fIn(opp(i),bbRegion);
end

% STREAMING STEP
for i=1:9
fIn(i,:,: ) = circshift(fOut(i,:,: ), [0,cx(i),cy(i)]);
end

% MICROSCOPIC BOUNDARY CONDITIONS: West boundary condition (Zou/He BC)
fIn(2,1,channel) = fIn(4,1,channel) + 2/3rho(:,1,channel).ux(:,1,channel);
fIn(6,1,channel) = fIn(8,1,channel) + 1/2
(fIn(5,1,channel)-fIn(3,1,channel))+ …
1/6
rho(:,1,channel).ux(:,1,channel) + 1/2rho(:,1,channel).uy(:,1,channel);
fIn(9,1,channel) = fIn(7,1,channel) + 1/2
(fIn(3,1,channel)-fIn(5,1,channel))+ …
1/6*rho(:,1,channel).ux(:,1,channel) - 1/2rho(:,1,channel).*uy(:,1,channel);

% % MICROSCOPIC BOUNDARY CONDITIONS: East boundary condition
fIn(7,lx,channel)=2fIn(7,lx-1,channel)-fIn(7,lx-2,channel);
fIn(4,lx,channel)=2
fIn(4,lx-1,channel)-fIn(4,lx-2,channel);
fIn(8,lx,channel)=2*fIn(8,lx-1,channel)-fIn(8,lx-2,channel);

% VISUALIZATION
if (mod(cycle,tPlot)==0)
u = reshape(sqrt(ux.^2+uy.^2),lx,ly);
u(bbRegion) = nan;
imagesc(u(:,ly:-1:1)’./uChannel);
colorbar
axis equal off; drawnow
end

end
toc;