Hi,
lobnat, I am also trying to understand half-way
and full-way bounce-back and I feel I’ve finally
grasped it. I commented out the Zou/He full-way BB
in the cavity.m and implemented the half-way
moving BB according to Ladd (1994) and Aidun
(1998). According to Ladd (2002), it is also
possible to replace the fluid local density with
the fluid mean density (i.e. unity).
Please have a look at the MATLAB code posted
below. I have run it and the results look fine.
Interested LB forum users, please let me know if I
made any mistake and I really wish more
discussions about half-way and full-way
bounce-back rules especially regarding moving
boundaries and (un)steady state flow (e.g.
particle suspensions).
BTW, in the code, East (West) Wall is actually
West (East) Wall and I wish I had named ewcol
wwrow, wwcol ewrow, swrow swcol and nwrow nwcol
instead. 
clear all
close all
clc
% GENERAL FLOW CONSTANTS
lx = 128;
ly = 128;
uLid = 0.05; % horizontal lid velocity
vLid = 0; % vertical lid velocity
Re = 100; % Reynolds number
nu = uLid lx / Re; % kinematic viscosity
omega = 1. / (3nu+1./2.); % relaxation parameter
maxT = 40000; % total number of iterations
tPlot = 10; % cycles for graphical output
% D2Q9 LATTICE CONSTANTS
% D2Q9 Lattice constants
% 7 3 6
% \ | / ^
% 4 - - 1 - - 2 | * * *
% / | \ | * * *
% 8 5 9 y x - ->
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];
lid = [2: (lx-1)];
= meshgrid(1:ly,1:lx);
obst = ones(lx,ly);
obst(2:(lx-1),2:(ly-1)) = 0; % Leave EWNS walls as
bbRegion
bbRegion = find(obst);
% 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,lxly)),
1,lx,ly ) ./rho;
uy = reshape ( (cy * reshape(fIn,9,lxly)),
1,lx,ly ) ./rho;
% % MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS
%
% ux(:,lid,ly) = uLid; %lid x - velocity
% uy(:,lid,ly) = vLid; %lid y - velocity
% rho(:,lid,ly) = 1 ./ (1+uy(:,lid,ly)) .* ( …
% sum(fIn([1,2,4],lid,ly)) +
2*sum(fIn([3,6,7],lid,ly)) );
% % MICROSCOPIC BOUNDARY CONDITIONS: LID (Zou/He
BC)
% fIn(5,lid,ly) = fIn(3,lid,ly) -
2/3rho(:,lid,ly).uy(:,lid,ly);
% fIn(9,lid,ly) = fIn(7,lid,ly) +
1/2(fIn(4,lid,ly)-fIn(2,lid,ly))+ …
% 1/2rho(:,lid,ly).*ux(:,lid,ly)
- 1/6rho(:,lid,ly).uy(:,lid,ly);
% fIn(8,lid,ly) = fIn(6,lid,ly) +
1/2(fIn(2,lid,ly)-fIn(4,lid,ly))- …
% 1/2rho(:,lid,ly).*ux(:,lid,ly)
- 1/6*rho(:,lid,ly).*uy(:,lid,ly);
% 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
% Mid-Wall Halfway Bounce-Back
% References:
% 1. Ladd, J. Fluid Mech. 1994a
% 2. Aidun et al., J. Fluid Mech. 1998
% East Wall x = 1, y = 2:(ly-1)
ewcol = 2:(ly-1);
fOut(2,1,ewcol) = fOut(4,2,ewcol);
fOut(6,1,ewcol) = fOut(8,2,ewcol+1);
fOut(9,1,ewcol) = fOut(7,2,ewcol-1);
% East Corner down x = 1, y = 1
fOut(6,1,1) = fOut(8,2,2);
% East Corner up x = 1, y = ly
rho(:,2,ly-1) =
fOut(1,2,ly-1)+fOut(3,2,ly-1)+fOut(4,2,ly-1)+fOut(
6,2,ly-1)…
+fOut(7,2,ly-1)+fOut(8,2,ly-1)+fOut(3,2,ly-1-1)+fO
ut(4,2+1,ly-1)+fOut(7,2+1,ly-1-1);
fOut(9,1,ly) = fOut(7,2,ly-1) +
2*rho(:,2,ly-1)3t(9)*uLid; % moving Bounce-Back
% West Wall x = lx, y = 2:(ly-1)
wwcol = 2:(ly-1);
fOut(4,lx,wwcol) = fOut(2,lx-1,wwcol);
fOut(8,lx,wwcol) = fOut(6,lx-1,wwcol-1);
fOut(7,lx,wwcol) = fOut(9,lx-1,wwcol+1);
% West Corner down x = lx, y = 1
fOut(7,lx,1) = fOut(9,lx-1,2);
% West Corner up x = lx, y = ly
rho(:,lx-1,ly-1) =
fOut(1,lx-1,ly-1)+fOut(2,lx-1,ly-1)+fOut(3,lx-1,ly
-1)…
+fOut(6,lx-1,ly-1)+fOut(7,lx-1,ly-1)+fOut(9,lx-1,l
y-1)…
+fOut(2,lx-1-1,ly-1)+fOut(3,lx-1,ly-1-1)+fOut(6,lx
-1-1,ly-1-1);
fOut(8,lx,ly) = fOut(6,lx-1,ly-1) +
2rho(:,lx-1,ly-1)3t(8)(-uLid); % moving
Bounce-Back
% South Wall x = 2:(lx-1), y =1
swrow = 2:(lx-1);
fOut(3,swrow,1) = fOut(5,swrow,2);
fOut(6,swrow,1) = fOut(8,swrow+1,2);
fOut(7,swrow,1) = fOut(9,swrow-1,2);
% North Wall x = 2:(lx-1), y = ly
nwrow = 2:(lx-1);
fOut(5,nwrow,ly) = fOut(3,nwrow,ly-1);
rho(:,3:(lx-2),ly-1) =
fOut(1,3:(lx-2),ly-1)+fOut(3,3:(lx-2),ly-1)…
+fOut(6,3:(lx-2),ly-1)+fOut(7,3:(lx-2),ly-1)+fOut(
2,(3:(lx-2))-1,ly-1)…
+fOut(4,(3:(lx-2))+1,ly-1)+fOut(3,3:(lx-2),ly-1-1)
+fOut(6,(3:(lx-2))-1,ly-1-1)…
+fOut(7,(3:(lx-2))+1,ly-1-1);
fOut(8,nwrow,ly) = fOut(6,nwrow-1,ly-1) +
2.*rho(:,nwrow-1,ly-1).*3.t(8).(-uLid); % moving
Bounce-Back
fOut(9,nwrow,ly) = fOut(7,nwrow+1,ly-1) +
2.*rho(:,nwrow+1,ly-1).*3.*t(9).*uLid; % moving
Bounce-Back
% % MICROSCOPIC BOUNDARY CONDITIONS: NO-SLIP WALLS
(bounce-back) Mid-wall
% Fullway 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)’./uLid);
colorbar
axis equal off; drawnow
end
end
Thanks you for the post.
Hi guys, Im a newbie. Nice to join this forum.