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,lx*ly)),*

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/3*rho(:,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/6
*rho(:,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)*3*t(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) +

2*rho(:,lx-1,ly-1)**3*t(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.