I have changed a single-relaxation-time matlab code into a two-relxation-time one. The final volocity contour seems abnormal. The Matlab code
clear
close all
clc
nx = 101;
ny = 101;
f = zeros(nx,ny,9);
feq = zeros(nx,ny,9);
u = zeros(nx,ny);
v = zeros(nx,ny);
rho = ones(nx,ny);
weight = [1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36 4/9];
cx = [1 0 -1 0 1 -1 -1 1 0];
cy = [0 1 0 -1 1 1 -1 -1 0];
opp = [3 4 1 2 7 8 5 6 9];
U = 0.020;
tau1 = 1.2;
omega1 = 1 / tau1;
tau2 = 0.25/(tau1-0.5)+0.5;
omega2 = 1/tau2;
w1 = 0.5*(omega1+omega2);
w2 = 0.5*(omega1-omega2);
count = 0;
tol = 1.0e-6;
error = 10.;
erso = 0.0;
%Main Loop
h = animatedline('Color','b','LineStyle','-','LineWidth',1);
set(gca,'YScale','log');
while error > tol
% Collitions
[f] = collision(nx,ny,u,v,cx,cy,weight,w1, w2,f,rho, opp);
% Streaming:
[f] = stream(f);
% End of streaming
%Boundary condition:
[f] = boundary(nx,ny,f,U,rho);
% Calculate rho, u, v
[rho,u,v] = ruv(nx,ny,f);
count = count+1;
ers = 0.;
for i = 1:nx
for j = 1:ny
ers = ers+u(i,j)*u(i,j)+v(i,j)*v(i,j);
end
end
error = abs(ers-erso);
erso = ers;
% Residual monitoring
if mod(count,100) == 0
addpoints(h,count,error);
hold on
drawnow
end
end
figure()
contourf(u')
% Collition
function [f] = collision(nx,ny,u,v,cx,cy,weight,w1,w2,f,rho,opp)
for j = 1:ny
for i = 1:nx
t1 = u(i,j)*u(i,j)+v(i,j)*v(i,j);
for k = 1:9
t2 = u(i,j)*cx(k)+v(i,j)*cy(k);
feq(i,j,k) = rho(i,j)*weight(k)*(1.0+3.0*t2+4.5*t2*t2-1.5*t1);
end
for k = 1:9
n = opp(k);
f(i,j,k) = f(i,j,k)-w1*( f(i,j,k)-feq(i,j,k) )-w2*( f(i,j,n)-feq(i,j,n) );
end
end
end
end
% Streaming:
function [f] = stream(f)
f(:,:,1) = circshift( squeeze(f(:,:,1)), [+1,+0] );
f(:,:,2) = circshift( squeeze(f(:,:,2)), [+0,+1] );
f(:,:,3) = circshift( squeeze(f(:,:,3)), [-1,+0] );
f(:,:,4) = circshift( squeeze(f(:,:,4)), [+0,-1] );
f(:,:,5) = circshift( squeeze(f(:,:,5)), [+1,+1] );
f(:,:,6) = circshift( squeeze(f(:,:,6)), [-1,+1] );
f(:,:,7) = circshift( squeeze(f(:,:,7)), [-1,-1] );
f(:,:,8) = circshift( squeeze(f(:,:,8)), [+1,-1] );
end
function [f] = boundary(nx,ny,f,U,rho)
%Boundary condition:
%left boundary, u
for j = 1:ny
f(1,j,1) = f(1,j,3)+2./3.*rho(1,j)*U;
f(1,j,5) = f(1,j,7)-0.5*(f(1,j,2)-f(1,j,4))+rho(1,j)*U/6;
f(1,j,8) = f(1,j,6)+0.5*(f(1,j,2)-f(1,j,4))+rho(1,j)*U/6;
u(1,j)= U;
end
%right hand boundary
f(nx,:,3) = f(nx-1,:,3);
f(nx,:,7) = f(nx-1,:,7);
f(nx,:,6) = f(nx-1,:,6);
%bottom boundary, bounce back
f(:,1,2) = f(:,1,4);
f(:,1,5) = f(:,1,7);
f(:,1,6) = f(:,1,8);
%Top boundary, bounce back
f(:,ny,4) = f(:,ny,2);
f(:,ny,8) = f(:,ny,6);
f(:,ny,7) = f(:,ny,5);
end
function[rho,u,v] = ruv(nx,ny,f)
rho = sum (f,3);
for i = 1:nx
rho(i,ny) = f(i,ny,9)+f(i,ny,1)+f(i,ny,3)+2.*(f(i,ny,2)+f(i,ny,6)+f(i,ny,5));
end
%calculate velocity compnents
u = ( sum(f(:,:,[1 5 8]),3) - sum(f(:,:,[3 6 7]),3) )./rho;
v = ( sum(f(:,:,[2 5 6]),3) - sum(f(:,:,[4 7 8]),3) )./rho;
end
The velocity contour
Is this result correct? How to fix the problem? Any suggestions will be thanked.