Hi all,

From a code (lid driven cavity) i found on this forum (sorry I don’t remember the name of the user), I try to implement MRT model.

The first test I made is a comparison between LBGK and MRT with s_i = 1/tau. I read, in a lot of talks and papers that is this case I must have exactly the same results. But this is not the case for long time integration. For example for long times, the velocity is constant with BGK model but envoles (slightly) with MRT model and results differ a lot. I don’t think there is a mistake on the implementation, so I wonder if there is a special condition to obtain same results with MRT(s_i = 1/tau) and LBGK models or if it’s always true ! (whatever the boundaries conditions )

Thanks

```
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 2D lid driven cavity by using LBE Multi-Relaxation time scheme %%%
%%% Writer : Arman Safdari %%%
%%% E-mail : sarman2@live.utm.my %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;close all;clc;
%%%% definitions %%%%%%%
lx=101;
ly=101;
rho(1:lx,1:ly)=1;
Uo=-0.05;
Re = 100;
Nu=abs(Uo)*lx/Re;
tau=3*Nu+0.5;
u=zeros(lx,ly);
v=zeros(lx,ly);
fEq=ones(9,lx,ly);
mEq=ones(9,lx,ly);
f=ones(9,lx,ly);
t = [4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36];
ex=[0,1,0,-1,0,1,-1,-1,1];
ey=[0,0,1,0,-1,1,1,-1,-1];
S=[1,1.4,1.4,1,1.2,1,1.2,1/tau,1/tau];
omega=1/tau;
S(:)=1/tau; % bgk case
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];
InvMS=M\diag(S);
fid = fopen('lid_driven_stat.dat','w');
fprintf(fid,'#iter - max(u) - max(v) - min(v) \n');
u(1:lx,1:ly)=0;v(1:lx,1:ly)=0;rho(1:lx,1:ly)=1;
for i=1:lx
for j=1:ly
f(1:9,i,j) = t(1:9);
end
end
for time=1:50000
%%
%%%%%%% collision %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
test=1;
if ( test == 1 )
% MRT case
for i=1:lx
for j=1:ly
jx=rho(i,j) * u(i,j);jy=rho(i,j) * v(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
m = reshape(M * reshape(f,9,lx*ly),9,lx,ly);
f = f - reshape(InvMS * reshape(m-mEq,9,lx*ly),9,lx,ly);
elseif ( test == 2 )
% BGK case
for i=1:9
cu = 3*(ex(i)*u+ey(i)*v);
fEq(i,:,:) = rho .* t(i) .* ...
( 1 + cu + 1/2*(cu.*cu) - 3/2*(u.^2+v.^2) );
end
f = f - omega .* (f-fEq);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%% streaming %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:9
f(i,:,:) = circshift(f(i,:,:), [0,ex(i),ey(i)]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%% Boundary Condition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Left Wall (bounce back)
f(2,1,:)=f(4,1,:);
f(6,1,:)=f(8,1,:);
f(9,1,:)=f(7,1,:);
%%%% Right Wall (bounce back)
f(4,lx,:)=f(2,lx,:);
f(8,lx,:)=f(6,lx,:);
f(7,lx,:)=f(9,lx,:);
%%%% Bottom Wall (bounce back)
f(3,:,1)=f(5,:,1);
f(6,:,1)=f(8,:,1);
f(7,:,1)=f(9,:,1);
%%%% Top Wall (Zou/He)
rhoLid=f(1,:,ly)+f(2,:,ly)+f(4,:,ly)+2*(f(3,:,ly)+f(7,:,ly)+f(6,:,ly));
f(5,:,ly)=f(3,:,ly);
f(9,:,ly)=f(7,:,ly) + 1/2*(f(4,:,ly)-f(2,:,ly)) + (1/2 * rhoLid * Uo);
f(8,:,ly)=f(6,:,ly) + 1/2*(f(2,:,ly)-f(4,:,ly)) - (1/2 * rhoLid * Uo);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%% rho and Velocity %%%%%%%%%%%%%%%%%%%%%%%%%%%
rho = reshape(sum(f),lx,ly);
u(:,:) = reshape(ex * reshape(f,9,lx*ly),lx,ly)./rho;
v(:,:) = reshape(ey * reshape(f,9,lx*ly),lx,ly)./rho;
u(1,:) = 0 ; u(lx,:)= 0 ; u(:,1) = 0 ; u(:,ly) = Uo;
v(1,:) = 0 ; v(lx,:)= 0 ; v(:,1) = 0 ; v(:,ly) = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Plotting part %%%%%%%%%%%%%%%%%%%%%%%%
% Psi=zeros(lx,ly);
% for i =2:lx-1
% Psi(i,2:ly-1) = Psi(i-1,2:ly-1) - v(i,2:ly-1);
% end
% contour(0:1/(lx-1):1,0:1/(ly-1):1,Psi',25);
% drawnow
%disp(t)
if (mod(time,100) == 0 )
fprintf(fid,'%8.0f %12.8f %12.8f %12.8f\n',time,max(max(u(:,:)))/abs(Uo),max(max(v(:,:)))/abs(Uo),min(min(v(:,:)))/abs(Uo));
end
end
%%
%%%%%%%%%%%%%%% Plotting part %%%%%%%%%%%%%%%%%%%%%%%%
Psi=zeros(lx,ly);
for i =2:lx-1
Psi(i,2:ly-1) = Psi(i-1,2:ly-1) - v(i,2:ly-1);
end
contour(0:1/(lx-1):1,0:1/(ly-1):1,Psi',25);
drawnow
```