recovery of LBGK with MRT model

Lattice Boltzmann Method LBM: Theory
#1

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
#2

Hi,

MRT should most certainly reduce to BGK if all the relaxation times are set to be the same. We can see this quite easily - just multiply the matrices yourself. This means that
something is doing something it shouldn’t in your code. I haven’t looked through it but if you can’t find any other bugs, my guess would be it has something to do with conservation
trouble. If you set the density and momentum moments to have a (non-zero) equilibrium then of course there should be no “relaxation” because they are conserved (their equilibrium are themselves). However, if you are loosing mass or momentum somewhere (as can often be the case in lid driven cavity because there are stress singularities in the corner, and because boundary conditions can be troublesome) then there will be an (unphysical) relaxation. It may be better to set the conserved equilibria to be zero (and perhaps to monitor any mass or momentum loss in the system, and the accuracy, with BGK). That’s just a thought, but we can certainly show the equivalence.

Good luck!

#3

Hi pleb01,

thank for your answer ! I thought about it and my conclusion is that the problem comes form the collision term modelling. The corner singularities are the same for both BGK and MRT models so, in my opinion, it can not be the key!

I thing the explanation comes form the MRT model. As it’s described in Lallemand and Luo paper, the “density fluctuation is decoupled from the momentum equation, similar to an incompressible LBE model”. The associated equilibrium function is then


f_i^eq = w_i ( rho + rho0 ( 3 (c_i u) + 4.5 (c_i u)[sup]2[/sup] - 1.5 u[sup]2[/sup]  ) )

where rho0 is usually set to be 1.

My guess is that with the MRT model you have to use rho = rho0 in the calculation of u and v and obviously if zou/he formulation is used for boundary conditions you have to take care about it.
I’m not still sure about that because it is not written explicitly in the paper and in most articles i read people using MRT write


rho u = sum c_i f_i

without stipulate if rho = rho0 = 1 like in an incompressible LBE model

Best regards,

Ben

#4

Hi Ben,

You should certainly make sure that you are comparing like-with-like, ie comparative equilibria. Again, this is easy to check: either do some hand calculations or use the equilibrium you quote above in your BGK code (you may want to be careful if you have Zou and He boundaries but bounce back shouldn’t be an issue - and if you are using Zou and He for the moving wall, you could instead just use an equilibrium boundary condition, as in Hou et al 1994 - which is the “benchmark” LB paper for this flow).

A quick point to note: when the flow is unsteady, the “incompressible” LB and the “standard” LB have a compressibility error of the same order.

Also, although I haven’t really looked through your code, I did notice you explicitly set the velocity on the walls “by hand”. I wouldn’t do this. The conditions should be satisfied by the boundary method being used, and if they are not then there is a failure with the method (for example, we know bounce back doesn’t give exactly zero wall velocity). Imposing things this way is arguably inconsistent, can deteriorate stability, and (perhaps relevant to your current problem) could mean that the momentum calculated via an equilibrium function is not necessarily the same as that imposed/forced “by hand”.

Tim

#5

Hi,

I changed the way boundary condition are taken care. I implement a proper version with the incompressible model and I finally recover that MRT with all coefficients equal 1/tau exactly leads to LBGK.
MRT with the “classical” parameters defined in Lallemand and Luo paper describe an incompressible fluid.

Thanks for this helpful discussion !

Ben