# LBM Algorithm

Hello,
I am Emad. I started to learn LBM. I have written some LBM codes in MATLAB, but some of them do not work!
Is there anybody who can help me by telling about LBM algorithm for programming?!
Best regards,

I wish you could come up with more specific questions. I will gladly try to understand what’s you’re problem if you give give us more clues.

Good Luck

Albert P

Dear Albert

Thank you very much.
Actually, I need an algorithm for programming LBM. Because I have wrote a MATLAB code but it does not work.
I think my algorithm for programming is wrong.
I need the algorithm like this:

1. defining parameters
2. LBM loop
2-1) define rho, v, etc
2-2) Streaming

I would be so sincere if you could help me.

if you are looking for generic information on how to do a code you should look up this page

Here you can find many important and fundamental information and also some matlab/c++ codes of some important features as SRT/MRT, binary liquids, porous media, and other advanced stuff like IBM.

I am sorry but I do not work with matlab (yet) so I won’t be able to help you with programming problems.

again, good luck.

Albert P

Thank you dear Albert…

I am new to LBM_MRT, and just made a code for Poiseulli flow. it is producing results but with some error. Also, it is working with simple bounce back but when i am applying second order bounce back, it is not producing any result but zero velocity. can any body suggest me any correction. thanks. Here is code;
clc; clear all; close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation of 2-Dimenional Poiselli Flow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic;
% Input parameters
Ly=11;
Lx=3;
tau1=1;
g=0.001;
% F=tau*g;
% Geometry
solid=zeros(Ly,Lx);
for j=1:Lx;
solid(1,j)=1;
solid(Ly,j)=1;

end
display(‘solid’)
solid;

% Initial Conditions

cx=[0.0,1.0,0.0,-1.0,0.0,1.0,-1.0,-1.0,1.0];
cy=[0.0,0.0,1.0,0.0,-1.0,1.0,1.0,-1.0,-1.0];
w=[4./9.,1./9.,1./9.,1./9.,1./9.,1./36.,1./36.,1./36.,1./36.];
u=zeros(Ly,Lx); v=zeros(Ly,Lx);
rho=ones(Ly,Lx); F=zeros(9,Ly,Lx);

for k=1:9;
f(k,:,:)=w(k)rho;
end
f;
for k=1:9;
F(k,:,:)=3
w(k)*cx(k)gtau1;
end
% Define ‘M’ matrix

M = [1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0;
-4.0 -1.0 -1.0 -1.0 -1.0 2.0 2.0 2.0 2.0;
4.0 -2.0 -2.0 -2.0 -2.0 1.0 1.0 1.0 1.0;
0.0 1.0 0.0 -1.0 0.0 1.0 -1.0 -1.0 1.0;
0.0 -2.0 0.0 2.0 0.0 1.0 -1.0 -1.0 1.0;
0.0 0.0 1.0 0.0 -1.0 1.0 1.0 -1.0 -1.0;
0.0 0.0 -2.0 0.0 2.0 1.0 1.0 -1.0 -1.0;
0.0 1.0 -1.0 1.0 -1.0 0.0 0.0 0.0 0.0;
0.0 0.0 0.0 0.0 0.0 1.0 -1.0 1.0 -1.0];

% ‘S’ matrix
S =[1.0,1.4,1.4,1.0,1.2,1.0,1.2,tau1,tau1];

RTM=inv(M)*diag(S); % Relaxation time for N-S Equation

% Main Loop

for ts=1:3000;
disp(ts)
rho = reshape(sum(f),Ly,Lx);
u(:, = reshape(cx * reshape(f,9,LxLy),Ly,Lx)./rho;
v(:, = reshape(cy * reshape(f,9,Lx
Ly),Ly,Lx)./rho;

``````% Calculate Equilibrium Function
for i=1:Ly;
for j=1:Lx;

u(i,j)=u(i,j)+g*tau1;
v(i,j)=v(i,j);
``````

%
jx=u(i,j)*rho(i,j);
jy=v(i,j)*rho(i,j);

``````        meq(1,i,j)= rho(i,j);
meq(2,i,j)= -2*rho(i,j)+3.*(jx.^2+jy.^2);
meq(3,i,j)= rho(i,j)-3.*(jx.^2+jy.^2);
meq(4,i,j)= jx;
meq(5,i,j)= -jx;
meq(6,i,j)= jy;
meq(7,i,j)= -jy;
meq(8,i,j)= jx.^2+jy.^2;
meq(9,i,j)= jx*jy;

end
end
``````

% Collision in Moment Space
% if (i==1 || i==Ly) or if (solid(i,j)
% % % Bounce Back
% temp=f(2,i,j); f(2,i,j)=f(4,i,j); f(4,i,j)=temp;
% temp=f(3,i,j); f(3,i,j)=f(5,i,j); f(5,i,j)=temp;
% temp=f(6,i,j); f(6,i,j)=f(8,i,j); f(8,i,j)=temp;
% temp=f(7,i,j); f(7,i,j)=f(9,i,j); f(9,i,j)=temp;
%
% else

`````` m=reshape(M*reshape(f,9,Lx*Ly),9,Ly,Lx);
f=f-reshape(RTM*(reshape((m-meq),9,Ly*Lx)),9,Ly,Lx);
``````

% end
% Streaming Function
for i=1:9
f(i,:, = circshift(f(i,:,:), [0,cy(i),cx(i)]);
end

``````    %%%% Bottom Wall (bounce back)
f(3,1,:)=f(5,1,:);
f(6,1,:)=f(8,1,:);
f(7,1,:)=f(9,1,:);
%%%% Top Wall (bounce back)
f(5,Ly,:)=f(3,Ly,:);
f(8,Ly,:)=f(6,Ly,:);
f(9,Ly,:)=f(7,Ly,:);
``````

% This the end of N-S Equation Solution
% Calculate Microscopic density rho, Velocity u, v
for i=1:Ly;
for j=1:Lx;
% rho(i,j)=0.0;
% u(i,j)=0.0;
% v(i,j)=0.0;
% % if~solid(i,j);
% for k=1:9;
% rho(i,j)=rho(i,j)+f(k,i,j);
% u(i,j)=u(i,j)+cx(k)*f(k,i,j);
% v(i,j)=v(i,j)+cy(k)*f(k,i,j);
% end
%
% u(i,j)=u(i,j)/rho(i,j);
% v(i,j)=v(i,j)/rho(i,j);
% end
x(i,j)=j;
y(i,j)=i;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

figure
quiver(x,y,u,v)
width=Ly-2;
figure
plot(y(:,1)-1.5-width/2,u(:,1),‘ro’)
hold on

% Poiseuille velocity profile as blue line
nu=1/3*(tau1-1/2);
plot(y(:,1)-1.5-width/2,g/(2nu)((width/2)^2-(y(:,1)-1.5-width/2).^2));
% axis ([-6 6 0 0.1]);
% axis auto
grid on;
% End of Program

toc