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,

Emad

Dear Emad,

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:

- defining parameters
- LBM loop

2-1) define rho, v, etc

2-2) Streaming

… - …

I would be so sincere if you could help me.

Emad,

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)

*g*tau1;

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,Lx*Ly),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/(2*nu)*((width/2)^2-(y(:,1)-1.5-width/2).^2));

% axis ([-6 6 0 0.1]);

% axis auto

grid on;

% End of Program

toc