I am trying to implement Two-Relaxation Time (TRT) model to compute gaseous flows in a long microchannel. I am following the paper “Two Relaxation time Lattice Boltzmann Model for rarefied gaseous flows” by J.A. Esfahani and A. Norouzi Link to the paper

The code seems to work fine for a low Knudsen number (0.0194). However, it seems to diverge for 0.194 and above, as tau_a (anti-symmetrical relaxation time) tends to diverge. Can anyone mention if there is anything wrong in my TRT algorithm implemented below or in the microchannel code? The grid arrangement and formulations have been taken from the aforementioned paper itself.

The MATLAB code is given here:

[code=“cpp”]

format long

n=210; %grid size along x-direction

m=21; %grid size along y-direction

iter=0;

rho=2*ones(n+2,m);
cx=[0 1 0 -1 0 1 -1 -1 1]; %Lattice Directions
cy=[0 0 1 0 -1 1 1 -1 -1];
u=zeros(n+2,m);
v=zeros(n+2,m);
f=zeros(9,n+2,m);
u_max=0.10; %maximum velocity
kn=0.0194;
tau_s=0.5+sqrt(6/pi)*m)/alph; %reynolds number

*kn*m; alph=(tau_s-0.5)/3; %dynamic viscosity re=(u_max

w=[4/9;1/9;1/9;1/9;1/9;1/36;1/36;1/36;1/36];

x=(1:n+2); %lattice locations - on grid nodes

y=(1:m)-0.5; %mid-grid

rhoo=1.0; %density at outlet

rhoi=1.4; %density at inlet

rho(2,:)=rhoi;

rho(n+1,:)=rhoo;

err_max=1.0;

opp=[1,4,5,2,3,8,9,6,7];

tau_a=(8

*tau_s-1)/(8*(2

*tau_s-1)); %anti-symmetric relaxation time*

f=f_init(u,v,cx,cy,w,n,m,rho); %Initialize feq to f

u_t=zeros(n+2,1);

u_b=zeros(n+2,1);

u_top=zeros(n+2,1);

u_bot=zeros(n+2,1);

rhoi_new=ones(m,1);

rhoo_new=ones(m,1);

f_symm=zeros(5,n+2,m);

f_asymm=zeros(5,n+2,m);

vel=[1 2 6 9 3];

r=3alph/(kn

f=f_init(u,v,cx,cy,w,n,m,rho); %Initialize feq to f

u_t=zeros(n+2,1);

u_b=zeros(n+2,1);

u_top=zeros(n+2,1);

u_bot=zeros(n+2,1);

rhoi_new=ones(m,1);

rhoo_new=ones(m,1);

f_symm=zeros(5,n+2,m);

f_asymm=zeros(5,n+2,m);

vel=[1 2 6 9 3];

r=3

*m*rhoo+3

*alph);*

while iter<4e4

iter=iter+1;

if rem(iter,100)==0

iter

err_max

% save(‘re500_258x258_antipar_ar04’,’-v6’);

end

feq=f_init(u,v,cx,cy,w,n,m,rho); %Calculate equilibrium dist function

f_symm(vel,:,:)=f(vel,:,:)+f(opp(vel),:,:)-feq(vel,:,:)-feq(opp(vel),:,:);

f_asymm(vel,:,:)=f(vel,:,:)-f(opp(vel),:,:)-feq(vel,:,:)+feq(opp(vel),:,:);

f(vel,:,:)=f(vel,:,:)-0.5f_symm(vel,:,:)/tau_s-0.5

while iter<4e4

iter=iter+1;

if rem(iter,100)==0

iter

err_max

% save(‘re500_258x258_antipar_ar04’,’-v6’);

end

feq=f_init(u,v,cx,cy,w,n,m,rho); %Calculate equilibrium dist function

f_symm(vel,:,:)=f(vel,:,:)+f(opp(vel),:,:)-feq(vel,:,:)-feq(opp(vel),:,:);

f_asymm(vel,:,:)=f(vel,:,:)-f(opp(vel),:,:)-feq(vel,:,:)+feq(opp(vel),:,:);

f(vel,:,:)=f(vel,:,:)-0.5

*f_asymm(vel,:,:)/tau_a; %Collision*

f(opp(vel),:,:)=f(opp(vel),:,:)-0.5f_symm(vel,:,:)/tau_s+0.5

f(opp(vel),:,:)=f(opp(vel),:,:)-0.5

*f_asymm(vel,:,:)/tau_a;*

%Extrapolation inlet/outlet boundary conditions

f(vel(2:4),1,:)=2f(vel(2:4),2,:)-f(vel(2:4),3,:); %extrapolated inlet bc

%Extrapolation inlet/outlet boundary conditions

f(vel(2:4),1,:)=2

f(opp(vel(2:4)),n+2,:)=2

*f(opp(vel(2:4)),n-1,:)-f(opp(vel(2:4)),n-2,:); %extrapolated outlet bc*

% rho(2,:)=2rho(3,:)-rho(4,:); %inlet density

% rho(2,:)=2

% rho(n+1,:)=2

*rho(n-1,:)-rho(n-2,:); %outlet density*

sum_rhoi=sum(rho(2,:),2); %averaging densities at inlet/outlet

sum_rhoo=sum(rho(n+1,:),2);

rhoi_new(:)=rho(2,

rhoo_new(:)=rho(n+1,

f(2,2:n+2,:)=f(2,1:n+1,:); %Direction (1,0) - Streaming

f(3,:,2:m)=f(3,:,1:m-1); %Direction (0,1)

f(4,1:n+1,:)=f(4,2:n+2,:); %Direction (-1,0)

f(5,:,1:m-1)=f(5,:,2:m); %Direction (0,-1)

f(6,2:n+2,2:m)=f(6,1:n+1,1:m-1); %Direction (1,1)

f(7,1:n+1,2:m)=f(7,2:n+2,1:m-1); %Direction (-1,1)

f(8,1:n+1,1:m-1)=f(8,2:n+2,2:m); %Direction (-1,-1)

f(9,2:n+2,1:m-1)=f(9,1:n+1,2:m); %Direction (1,-1)

%BSR boundary conditions at the top and bottom walls

f(8,:,m)=rf(6,:,m)+(1-r)

sum_rhoi=sum(rho(2,:),2); %averaging densities at inlet/outlet

sum_rhoo=sum(rho(n+1,:),2);

rhoi_new(:)=rho(2,

*m*rhoi/sum_rhoi;rhoo_new(:)=rho(n+1,

*m*rhoo/sum_rhoo;f(2,2:n+2,:)=f(2,1:n+1,:); %Direction (1,0) - Streaming

f(3,:,2:m)=f(3,:,1:m-1); %Direction (0,1)

f(4,1:n+1,:)=f(4,2:n+2,:); %Direction (-1,0)

f(5,:,1:m-1)=f(5,:,2:m); %Direction (0,-1)

f(6,2:n+2,2:m)=f(6,1:n+1,1:m-1); %Direction (1,1)

f(7,1:n+1,2:m)=f(7,2:n+2,1:m-1); %Direction (-1,1)

f(8,1:n+1,1:m-1)=f(8,2:n+2,2:m); %Direction (-1,-1)

f(9,2:n+2,1:m-1)=f(9,1:n+1,2:m); %Direction (1,-1)

%BSR boundary conditions at the top and bottom walls

f(8,:,m)=r

*f(7,:,m); %Top wall*

f(5,:,m)=f(3,:,m);

f(9,:,m)=rf(7,:,m)+(1-r)

f(5,:,m)=f(3,:,m);

f(9,:,m)=r

*f(6,:,m);*

f(7,:,1)=rf(9,:,1)+(1-r)

f(7,:,1)=r

*f(8,:,1); %bottom wall*

f(3,:,1)=f(5,:,1);

f(6,:,1)=rf(8,:,1)+(1-r)

f(3,:,1)=f(5,:,1);

f(6,:,1)=r

*f(9,:,1);*

rho=reshape(sum(f),n+2,m); %compute fields

rho(2,:)=rhoi_new(:);

rho(n+1,:)=rhoo_new(:);

un=reshape(sum(f([2,6,9],:,:))-sum(f([4,7,8],:,:)),n+2,m)./rho;

vn=reshape(sum(f([6,3,7],:,:))-sum(f([8,5,9],:,:)),n+2,m)./rho;

u_top(:,1)=1.5u(:,m)-0.5

rho=reshape(sum(f),n+2,m); %compute fields

rho(2,:)=rhoi_new(:);

rho(n+1,:)=rhoo_new(:);

un=reshape(sum(f([2,6,9],:,:))-sum(f([4,7,8],:,:)),n+2,m)./rho;

vn=reshape(sum(f([6,3,7],:,:))-sum(f([8,5,9],:,:)),n+2,m)./rho;

u_top(:,1)=1.5

*u(:,m-1);*

u_bot(:,1)=1.5u(:,1)-0.5*u(:,2);

u_bot(:,1)=1.5

t1=reshape(un-u,[],1); t2=reshape(vn-v,[],1);

A=u.^2+v.^2;

err_max=sqrt(t1’

*t1+t2’*(u_bot-u_b))/sqrt(sum(A(:))+sum(u_bot(:))+sum(u_top(:))); %computing residual error

*t2+(u_top-u_t)’*(u_top-u_t)+(u_bot-u_b)’u=un;v=vn; u_t=u_top; u_b=u_bot;

end

u_c=max(abs(u(n,:)));

u_new=cat(2,u_bot(:),u);

u_new=[u_new u_top(:)];

y_new=zeros(1,m+2);

y_new(2:m+1)=y(:);

y_new(m+2)=m;

plot(y_new/m,u_new(n,:)/u_c);

```
[code="cpp"]
function [ f ] = f_init(u,v,cx,cy,w,n,m,rho)
%initialisation of f to feq at t=0
f=zeros(9,n+2,m);
t1=u.^2+v.^2; %sum of squares of velocities
for k=1:9
t2=u*cx(k)+v*cy(k);
f(k,:,:)=rho.*w(k).*(1.0+3.0*t2+4.5*(t2.*t2)-1.50*t1);
end
end
```

Thanks

Sthavishtha