Moving Boundary using Half way bounce back!

Hello all,

I am a new beginner in the Lattice Boltzmann community and trying to learn my way through by simulating classical validation cases of Poiseuilli flow and lid driven cavity. I am able to get satisfactorily results for Poiseuilli flow however I am not sure how to use half way bounce back approach for the moving boundaries (lid driven cavity).

I am trying to follow: Momentum transfer of a Boltzmann-Lattice fluid with boundaries", Lallemand et. al., 2001. But I am not able to implement it properly as the results are not coherent with the benchmarks.

I am not sure but I think the problem lies in calculation of the gradient in eq 8 of the above paper. I was wondering if someone can shed some light on the aspect!

Also, do we need to use the equilibrium distribution function mentioned in an eq 10 and eq 11.

Hello,

first off, have you successfully implemented fullway bounceback for the resting and moving walls? I think it is a bit easier for the beginning. After you fully understand it, you can switch to halfway bounceback.

Timm

I have validated the full and half way bounce back for resting walls. But I haven’t got any success with the moving boundaries.

Any suggestions for the moving boundaries and full bounce back (other than the Zou/He bc).

Thanks again!

Hello again,

I try to explain the fullway bounceback and how it may be implemented (I have never done halfway bounceback).

This is a sketch of an upper wall in 2D (F = fluid, W = wall):


W W W W W
\|/
F F O F F



Populations are sent from the middle fluid node “O” to the three neighboring wall nodes. This happens at time step t. Let the populations be called 3 (top), 5 (top right), and 7 (top left). In the fullway bounceback, those populations remain on the wall until the next time step at t+1. At t+1, they are sent back to the original fluid node (opposite direction). Let those populations be called 3’, 5’, and 7’. Note that there is no collision on wall nodes. Collision is replaced by bounceback, i.e., by swapping the lattice velocities.
Now, one can show that, if the wall moves with velocity v to the right (or the left), the populations 3’, 5’, and 7’ at time t+1 are related to the populations 3, 5, and 7 at time t like this:
$f’_i = f_i - 2 \rho w_i (\vec c_i \cdot \vec v)/ c_s^2$
where w_i are the lattice weights, \vec c_i are the lattice velocities, c_s is the speed of sound and \vec v is the velocity vector of the wall. \rho is the density at fluid node “O” at time t.

This is all to correctly implement moving wall bounceback. It is exactly mass conserving, and it is second-order accurate if the wall is located in the middle between the fluid and wall nodes.

Timm

Many thanks Timm for the detailed description!

I seek a confirmation in your explanation. In addition to the special treatment of the directional density component [3’ 5’ 7’] we also need to swap the remaining components without addition/substraction of the momentum vector. As in full bounce back we need to reverse all the components!
i.e.
f’_i(1)=f_i(1);
f’_i(2)=f_i(2);
f’_i(4)=f_i(4); and so on…

Thanks in anticipation!

There is no need to bounceback the other populations. Why? Because they would be reflected into the wall where you do not have any fluid. Only populations going back into the fluid have to be taken care of. Populations trapped in the wall will never reappear, at least as long as fluid and wall nodes are fixed during the simulation.

Hello Tim,

Many thanks for your explanation again!

But as per my understanding, the reversal of all 9 components should not affect the solution. I am clueless about the problem.

I am certain about your busy schedule but if your time permits, may I send you my code, written in matlab, and if you can have a look at it!

Thanks in anticipation!

Maybe you should post the bounceback part of your code here so that everybody can have a look at it.

Hello Everyone!

A quick recapitulation of the discussion till now:-

1. I am facing problems with moving wall conditions (Bounce back!) Though I understood well the Zou/He condition and was able to validate it.
2. Timm suggested me to use full bounce back with addition/substraction of momentum in directional density components!
3. Using Timm’s suggestions, the effect of moving wall on flow field is observed. The code is stable but it is not giving the desired velocity profiles!

I could not figure out the possible cause of error and wondering if some of you can kindly devote some time and have a look at the implementation of the wall boundary condition. Look forward to hear some inputs from this community!

Thanks and regards!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X increasing from left to right
% Y increasing from top to bottom
% Top wall moving towards right
% D2Q9 Lattice Structure
%7 3 6
% \ | /
%4-1-2
% / |
%8 5 9
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Time step begins

% CALCULATION OF HYDROVARIABLES
ux(:, = (fIn(:,:,2) - fIn(:,:,4) + fIn(:,:,6) - fIn(:,:,7) - fIn(:,:,8) + fIn(:,:,9))./rho(:,:); %ux = u velocity
uy(:, = (fIn(:,:,6) + fIn(:,:,3) + fIn(:,:,7)- fIn(:,:,8) - fIn(:,:,5) - fIn(:,:,9))./rho(:, ; %uy = v velocity

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

% WALL BOUNDARY
wallX=2:lx-1;
for i=1:9 % FULL BOUNCE BACK (Avoid Collision by using fIn and not ftemp)
ftemp(:,1 ,i) = fIn(:,1 ,opp(i)); % WEST
ftemp(:,lx,i) = fIn(:,lx,opp(i)); % EAST
ftemp(ly,:,i) = fIn(ly,:,opp(i)); % SOUTH
ftemp(1,wallX,i) = fIn(1,wallX,opp(i)); % NORTH (First Reversal)
end
% NORTH TIMM’S MOVING WALL (Addition/Substraction of Momentum
for i=[5 8 9]
ftemp(1,wallX,i) = ftemp(1,wallX,i)+2*t(i)*cx(i)*uLid/oneothree;
end
fIn=ftemp;

% STREAMING STEP
fIn( : , 2:lx , 2) = ftemp( : , 1:lx-1, 2);
fIn(1:ly-1, : , 3) = ftemp(2:ly , : , 3);
fIn( : , 1:lx-1, 4) = ftemp( : , 2:lx , 4);
fIn(2:ly , : , 5) = ftemp(1:ly-1, : , 5);
fIn(1:ly-1, 2:lx , 6) = ftemp(2:ly , 1:lx-1, 6);
fIn(1:ly-1, 1:lx-1, 7) = ftemp(2:ly , 2:lx , 7);
fIn(2:ly , 1:lx-1, 8) = ftemp(1:ly-1, 2:lx , 8);
fIn(2:ly , 2:lx , 9) = ftemp(1:ly-1, 1:lx-1, 9);

% Time step ends here

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. / (3*nu+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)];

[y,x] = 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,lx*ly)), 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/2*rho(:,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/2*rho(:,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)+fOut(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,ly-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



Hello Jonny!

Many thanks for sharing your knowledge base and code.

I am trying to identify the difference between our implementation and will update you about my progress.

Thanks and regards!

Hi there
I’m also trying to move an object in a fluid.

I’ve implemented mid-plane full bounce back BB boundaries successfully (according to Sokup book), my case is very simple, and there is no curved boundaries. I’m trying to move a solid polygon inside the fluid.
anyway I dont know where to calculate $f’_i = f_i - 2 \rho w_i (\vec c_i \cdot \vec v)/ c_s^2$ that Timm suggested and if it would work for a solid moving in an arbitrary direction

also I have found on different documentation that ?rst-order equilibrium distribution at moving boundary can be calculated this way:

f^{(eq)}{i}i =f0{i} +?i * ei ·ew

ew: velocity of boundary
weight coef?cients

f0_{i} = {4/9,1/9,1/9,1/9,1/9, 1 /36, 1 /36, 1 /36, 1 /3
?i = {0,1/3,1/3,1/3,1/3, 1/ 12, 1/ 12, 1/ 12, 1/ 12}

but I’m not really sure about where to calculate this.

Hi there

Timm may I ask how did you end up with $f’_i = f_i - 2 \rho w_i (\vec c_i \cdot \vec v)/ c_s^2$??

according to the article (Lattice Boltzmann method for moving boundaries Pierre Lallemand, Li-Shi Luo)

after interpolation and substitution of q=1/2 (this is always going to happen for mid-plane bounce back) you get

f’_i = f_i + 3+w_a*(e_a + u_w)

where w_a are the weights and u_w is the speed of the wall

What is q and why interpolation? You are talking about interpolated bounceback, right? I am talking about bounceback for a plain wall without interpolation. The equation I have posted is correct. It is given in Ladd’s papers, and I have double-checked it.

Hello Jmazo,

I certainly can help you with the mid plane bounce back. The algorithm in nutshell is:

1. Complete the collision step on all nodes. (Use the general equilibrium distribution function. No need for a modified distribution function).
2. Do the streaming on all the nodes. (There are no ghost nodes so the wall is practically out of the domain).
3. Replace the unknown populations at the boundary fluid nodes with known populations in the opposite directions. Use the f values obtain after the collision.
4. Add/subtract the momentum vector in the bounce back population only (either the momentum value suggested by Timm or the one suggested in the Lallemand et. al). They are both the same for q=0.5).

@Timm: q is the distance between the fluid node and the wall.

I have validated it and it works perfectly fine.

Having said that, I am not sure how to move the polygon in the flow field. Do you calculate the net stress tensor acting on the polygon and then define the velocity of the polygon? It will be very interesting if you share your knowledge in this forum!

Thanks and regards!

I have found the formula $f’_i = f_i - 2 \rho w_i (\vec c_i \cdot \vec v)/ c_s^2$ in Ladd’s paper? but why it is?

You can for example derive the formula by transforming to a coordinate system where the wall is at rest (Galilei transform). In the resting wall system, f’_i = f_i holds (resting wall bounceback). You then have to transform back to the initial coordinate frame, and you automatically get the correct equation (up to second-order accuracy).

Hi lobnat,

In reference to your question involving moving polygons, and what I think jmazo is implementing, the paper by Lallemand and Luo, http://research.nianet.org/~luo/Reprints-luo/2003/Lallemand_JCP-2003.pdf, attempts moving objects within flows as well as ambient conditions.

jonny Wrote:

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
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

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. / (3
nu+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,lx
ly)),
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/2
rho(:,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/2
rho(:,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.

Watch Grown Ups Online Free