# Two Relaxation Time (TRT) Implementation for rarefied gaseous flows in microchannel

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=2ones(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)knm;
alph=(tau_s-0.5)/3; %dynamic viscosity
re=(u_max
m)/alph; %reynolds number
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=(8tau_s-1)/(8(2tau_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=3
alph/(knmrhoo+3alph);
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_symm(vel,:,:)/tau_s-0.5f_asymm(vel,:,:)/tau_a; %Collision
f(opp(vel),:,:)=f(opp(vel),:,:)-0.5
f_symm(vel,:,:)/tau_s+0.5f_asymm(vel,:,:)/tau_a;
%Extrapolation inlet/outlet boundary conditions
f(vel(2:4),1,:)=2
f(vel(2:4),2,:)-f(vel(2:4),3,:); %extrapolated inlet bc
f(opp(vel(2:4)),n+2,:)=2f(opp(vel(2:4)),n-1,:)-f(opp(vel(2:4)),n-2,:); %extrapolated outlet bc
% rho(2,:)=2
rho(3,:)-rho(4,:); %inlet density
% rho(n+1,:)=2rho(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,mrhoi/sum_rhoi;
rhoo_new(:)=rho(n+1,mrhoo/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(6,:,m)+(1-r)f(7,:,m); %Top wall
f(5,:,m)=f(3,:,m);
f(9,:,m)=r
f(7,:,m)+(1-r)f(6,:,m);
f(7,:,1)=r
f(9,:,1)+(1-r)f(8,:,1); %bottom wall
f(3,:,1)=f(5,:,1);
f(6,:,1)=r
f(8,:,1)+(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.5
u(:,m)-0.5u(:,m-1);
u_bot(:,1)=1.5
u(:,1)-0.5*u(:,2);
t1=reshape(un-u,[],1); t2=reshape(vn-v,[],1);
A=u.^2+v.^2;
err_max=sqrt(t1’t1+t2’t2+(u_top-u_t)’(u_top-u_t)+(u_bot-u_b)’(u_bot-u_b))/sqrt(sum(A(:))+sum(u_bot(:))+sum(u_top(:))); %computing residual error
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

Dear All

I managed to figure out the problem. Its pertaining to the boundary conditions, not the TRT implementation. However, I am facing some issues with the corner boundary conditions. The computational and physical domain of the channel is shown below The blue lines refer to the channel boundaries. https://docs.google.com/drawings/d/1q1HTElVT5TDGBxV7Kr5ZQwgwCeWPaTS1j67lZLB7HOI/pub?w=983&h=726

As per this figure, mid-grid nodes are used along the y-direction, and on-grid nodes along the x-direction. Lattice nodes outside the boundary along x-direction (1,3) are required for extrapolation (inlet) boundary conditions with fixed pressure. And the solid nodes along y-direction (3,5) are required for applying specular reflection boundary conditions (applied only for mid-grid nodes).

My doubt is which node should I consider as my corner node,where I need to apply the corner boundary conditions? What corner boundary conditions do I need to apply, as this is a mixture of mid-grid and on-grid nodes? Moreover, are the nodes (1,5) and (2,5) required/should they be omitted, if they don’t act as corner nodes?