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_maxm)/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=3alph/(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.5f_symm(vel,:,:)/tau_s-0.5f_asymm(vel,:,:)/tau_a; %Collision
f(opp(vel),:,:)=f(opp(vel),:,:)-0.5f_symm(vel,:,:)/tau_s+0.5f_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
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,:)=2rho(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)=rf(6,:,m)+(1-r)f(7,:,m); %Top wall
f(5,:,m)=f(3,:,m);
f(9,:,m)=rf(7,:,m)+(1-r)f(6,:,m);
f(7,:,1)=rf(9,:,1)+(1-r)f(8,:,1); %bottom wall
f(3,:,1)=f(5,:,1);
f(6,:,1)=rf(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.5u(:,m)-0.5u(:,m-1);
u_bot(:,1)=1.5u(:,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