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,:,:)=3w(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,LxLy),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