Welcome! Log In Create A New Profile

Advanced

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

Posted by sthavishthasthavishtha  
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:

Language: C++
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)*kn*m; 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=(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=3*alph/(kn*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.5*f_symm(vel,:,:)/tau_s-0.5*f_asymm(vel,:,:)/tau_a; %Collision f(opp(vel),:,:)=f(opp(vel),:,:)-0.5*f_symm(vel,:,:)/tau_s+0.5*f_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,:)=2*f(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,:)=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,:)*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(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.5*u(:,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);

Language: C++
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.

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?

I would be grateful for a quick reply. Thanks in advance.
Sorry, you do not have permission to post/reply in this forum.