I am trying to simulate a flow past cylinder at high Re using MRT Lbm. I am a newbie at this so i am not able to understand somethings. Firstly, calculating the equilibrium f values in normal lbm, we get such that RHO before is equal to Rho eqbm. So when applying to the formulat R-Req, wouldnt the density terms always cancel out?
That is correct … part of the definition for the local equilibrium distribution function is that its sum over all lattice links should equal the same as the sum of the distribution functions. No matter what the relaxation time for the density would be in MRT, the change in density due to collision is always zero and so the density will be the same both before and after the collision step.
The same is often true for the various momentum components, although the local equilibrium momentum may be adjusted by any body or interface force if these are applied. In this case, the relaxation time for momentum terms must be equal to 1 to ensure the correct fluid velocity is obtained.
i wrote a code to do the MRT-LBM around a cylinder but for some reason, it is not working. if you could help me figure out what is wrong, it would be greatly helpful.
the code is in two functions init() and main()
i think the problem is in the collision step with all the matrix manipulations, but im not sure yet.
the code is below.
function init()
global nx
global ny
global directions
global matrix
global ex
global ey
global tau
global tend
global M
global MI
global S
to = [4/9,1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
nx=100;
ny=50;%350;
density=1;
Uo=0.1;
Nu=0.01;
Re=abs(Uo)nx/Nu;
tau=3Nu+0.5;
M=[1, 1, 1, 1, 1, 1, 1, 1, 1;
-4,-1,-1,-1,-1, 2, 2, 2, 2;
4,-2,-2,-2,-2, 1, 1, 1, 1;
0, 1, 0,-1, 0, 1,-1,-1, 1;
0,-2, 0, 2, 0, 1,-1,-1, 1;
0, 0, 1, 0,-1, 1, 1,-1,-1;
0, 0,-2, 0, 2, 1, 1,-1,-1;
0, 1,-1, 1,-1, 0, 0, 0, 0;
0, 0, 0, 0, 0, 1,-1, 1,-1 ];
S=diag([1,1.4,1.4,1,1.2,1,1.2,1/tau,1/tau]);
MI=inv(M);
directions=9;
ex=[0,1,0,-1,0,1,-1,-1,1];
ey=[0,0,1,0,-1,1,1,-1,-1];
matrix=zeros(ny,nx,directions+1);
tend=1000;
obs_x=nx/5;
obs_y=ny/2;
obs_r=ny/6;
matrix([1 ny],:,10)=1;
for i=2:ny-1
for j=1:nx
if sqrt((i-obs_y)^2+(j-obs_x)^2)<=obs_r
matrix(i,j,10)=1;
end
for k=1:9
matrix(i,j,k)=density*to(k);
end
end
end
end
//////////////////////////////////////////////////////////////////////////////////
main function
/////////////////////////////////////////////////////////////////////////////////
clear all
close all
global nx
global ny
global directions
global matrix
global ex
global ey
global tau
global tend
global M
global MI
global S
tstart=tic;
init();
f1=3.0;
f2=9/2;
f3=3/2;
F=zeros(9,1);
R=zeros(9,1);
Req=R;
ux=zeros(ny,nx);
ueqx=zeros(ny,nx);
uy=zeros(ny,nx);
ueqy=zeros(ny,nx);
u=zeros(ny,nx);
feq=zeros(ny,nx,directions);
rho=zeros(ny,nx);
temp=zeros(directions+1,1);
MIS=MI*S;
for t=1:tend
disp(t);
%streaming step
for i=1:8
matrix(:,:,i)=circshift(matrix(:,:,i),[ey(i) ex(i)]);
end
%collision on walls (bounce back)
for i=1:ny
for j=1:nx
if matrix(i,j,10)==1 %if it is a wall
temp(1)=matrix(i,j,1);
temp(2)=matrix(i,j,4);
temp(3)=matrix(i,j,5);
temp(4)=matrix(i,j,2);
temp(5)=matrix(i,j,3);
temp(6)=matrix(i,j,8);
temp(7)=matrix(i,j,9);
temp(8)=matrix(i,j,6);
temp(9)=matrix(i,j,7);
temp(10)=1;
matrix(i,j,:)=temp(:);
end
end
end
% pressure boundary in right
for k=2:ny-1
fj=matrix(k,nx,:);
b=sum(fj);
ux0=-1+( fj(1)+fj(3)+fj(5)+2.*(fj(2)+fj(6)+fj(9)))/(b);
ru=(b)*ux0;
matrix(k,nx,4)=fj(2)-(2/3)*ru;
matrix(k,nx,8)=fj(6)-(1/6)ru+(1/2)(fj(3)-fj(5));
matrix(k,nx,7)=fj(9)-(1/6)ru+(1/2)(fj(5)-fj(3));
end
% velocity boundary in left
for k=2:ny-1
fj=matrix(k,1,:);
rho0=( fj(1)+fj(3)+fj(5)+2*(fj(4)+fj(8)+fj(7)))/(1-0.1);
ru=(0.1)*rho0;
matrix(k,1,2)=fj(4)+(2/3)*ru;
matrix(k,1,6)=fj(8)+(1/6)ru-(1/2)(fj(3)-fj(5));
matrix(k,1,9)=fj(7)+(1/6)ru-(1/2)(fj(5)-fj(3));
end
%calculating macro variables
rho=sum(matrix(:,:,:),3);
ux=(matrix(:,:,2)+matrix(:,:,6)+matrix(:,:,9)-matrix(:,:,4)-matrix(:,:,7)-matrix(:,:,8))./rho;
uy=(matrix(:,:,3)+matrix(:,:,6)+matrix(:,:,7)-matrix(:,:,5)-matrix(:,:,8)-matrix(:,:,9))./rho;
for i=1:ny
for j=1:nx
if matrix(i,j,10)==1
rho(i,j)=0;
ux(i,j)=0;
uy(i,j)=0;
end
end
end
%calculating equilibrium conditions
for i=1:ny
for j=1:nx
if matrix(i,j,10)==0
F=[matrix(i,j,1); matrix(i,j,2);matrix(i,j,3); matrix(i,j,4); matrix(i,j,5); matrix(i,j,6); matrix(i,j,7); matrix(i,j,8);matrix(i,j,9)];
rhoeq=rho(i,j);
jx=rho(i,j)*ux(i,j);
jy=rho(i,j)*uy(i,j);
eeq=-2*rho(i,j)+3*(jx^2+jy^2);
eeeq=rho(i,j)-3*(jx^2+jy^2);
jxeq=jx;
jyeq=jy;
qxeq=-1*jx;
qyeq=-1*jy;
pxxeq=jx^2-jy^2;
pxyeq=jx*jy;
Req=[rhoeq;eeq;eeeq;jxeq;qxeq;jyeq;qyeq;pxxeq;pxyeq];
R=M*F;
F=F-MIS*(R-Req);
matrix(i,j,1)=F(1);
matrix(i,j,2)=F(2);
matrix(i,j,3)=F(3);
matrix(i,j,4)=F(4);
matrix(i,j,5)=F(5);
matrix(i,j,6)=F(6);
matrix(i,j,7)=F(7);
matrix(i,j,8)=F(8);
matrix(i,j,9)=F(9);
end
end
end
%visualising
u=sqrt(ux.^2.+uy.^2);
imagesc(u);
axis equal off;
drawnow;
end
toc(tstart)
figure(2);
clf;
hold on;
axis equal off;
colormap(gray(2));
image(2-matrix(:,:,10));
quiver(ux,uy)
/////////////////////////////////////////////////////////////////////////////////////
thnx in advance