Hello Dear

Actually I have done my matlab code for Channel with Multi relaxation time but i think it has a problem with boundary condition and for example after 130000 iteration my code is crashed ,I was wondering If you could check my code and then say me what kind of boundary condition I should apply for Channel . I’m looking forward to hearing from you. Thank you for your consideration and time.

Yours faithfully

Mohammad Pourtousi

my matlab code:

clear;clc;

%%%% definitions %%%%%%%

lx=400;

ly=30;

hight1 =10;

hight2 =10;

lenght1=160;

lenght2=100;

lyInlet = hight1+1:ly;

lxInlet = 1:lenght1;

lxOutlet = lenght1+lenght2 : lx ;

lyOutlet = hight2 + 1 : ly ;

lyLeftCavity = 1:hight1;

lxBottomCavity = lenght1:lenght1+lenght2;

lyRightCavity = 1:hight2;

rho(1:lx+1,1:ly)=1;

Uo=0.05;

Re=100;

Nu = Uo*(ly-hight1)/Re;

tau = 3*Nu+0.5;

u=zeros(lx,ly);

v=zeros(lx,ly);

mEq=ones(9,lx,ly);

f=ones(9,lx,ly);

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

ex=[0,1,0,-1,0,1,-1,-1,1];

ey=[0,0,1,0,-1,1,1,-1,-1];

S=[1,1.4,1.4,1,1.2,1,1.2,1/tau,1/tau];

M=[1, 1, 1, 1, 1, 1, 1, 1, 1;-4,-1,-1,-1,-1, 2, 2, 2, 2;

4,-2,-2,-2,-2, 1, 1, 1, 1; 0, 1, 0,-1, 0, 1,-1,-1, 1;

0,-2, 0, 2, 0, 1,-1,-1, 1; 0, 0, 1, 0,-1, 1, 1,-1,-1;

0, 0,-2, 0, 2, 1, 1,-1,-1; 0, 1,-1, 1,-1, 0, 0, 0, 0;

0, 0, 0, 0, 0, 1,-1, 1,-1];

InvMS=M\diag(S);

for

t=1:200000

%%

%%%%%%% collision %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for

i=1:lx

for j=1:ly

jx=rho(i,j) * u(i,j);jy=rho(i,j) * v(i,j);

mEq(1,i,j)= rho(i,j);

mEq(2,i,j)= -2*rho(i,j) + 3 * (jx^2 + jy^2);

mEq(3,i,j)= rho(i,j) - 3 * (jx^2 + jy^2);

mEq(4,i,j)= jx;

mEq(5,i,j)= -jx;

mEq(6,i,j)= jy;

mEq(7,i,j)= -jy;

mEq(8,i,j)= jx^2 - jy^2;

mEq(9,i,j)= jx * jy;

end

end

m = reshape(M * reshape(f,9,lx*ly),9,lx,ly);

f = f - reshape(InvMS * reshape(m-mEq,9,lx*ly),9,lx,ly);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%

%%%%%%%% streaming %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for

i=1:9

f(i,:, = circshift(f(i,:,:), [0,ex(i),ey(i)]);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%

%%%%%%% Boundary Condition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% Left Wall (bounce back)

rhoInput=f(1,1,lyInlet)+f(3,1,lyInlet)+f(5,1,lyInlet)+2*(f(4,1,lyInlet)+

…

f(7,1,lyInlet)+f(8,1,lyInlet))/(1-Uo);

f(2,1,lyInlet)=f(4,1,lyInlet) + 2/3 * rhoInput * Uo;

f(6,1,lyInlet)=f(8,1,lyInlet) + rhoInput * Uo / 6 ;

f(9,1,lyInlet)=f(7,1,lyInlet) + rhoInput * Uo / 6 ;

% %%%% Right Wall (bounce back)

% f(2,1,lyInlet)= f(2,lx,lyOutlet) ;

% f(6,1,lyInlet)= f(6,lx,lyOutlet) ;

% f(9,1,lyInlet)= f(9,lx,lyOutlet);

%

% f(2,lx+1,lyOutlet)= f(2,1,lyInlet) ;

% f(6,lx+1,lyOutlet)= f(6,1,lyInlet) ;

% f(9,lx+1,lyOutlet)= f(9,1,lyInlet);

f(2,lx,lyOutlet)= f(2,1,lyInlet) ;

f(6,lx,lyOutlet)= f(6,1,lyInlet) ;

f(9,lx,lyOutlet)= f(9,1,lyInlet);

% -------------------------------------------------------------------------

% for i=1:9

% rho(lx,lyOutlet)= rho(lx-1,lyOutlet);

% v(lx,lyOutlet)=0;

% u(lx,lyOutlet)=u(lx-1,lyOutlet);

% % Right boundary

% % f(i,lx,lyOutlet) = mEq(i,lx,lyOutlet) + …

% % 18*w(i)*ex(i)*ey(i)* ( f(6,lx,lyOutlet) - …

% % f(9,lx,lyOutlet)-mEq(6,lx,lyOutlet)+mEq(9,lx,lyOutlet) );

% end

%%%% Bottom Wall (bounce back)

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

f(6,:,1)=f(8,:,1);

f(7,:,1)=f(9,:,1);

%%%% Top Wall

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

f(9,:,ly)=f(7,:,ly);

f(8,:,ly)=f(6,:,ly);

%%%% Bottom First Obstacle wall

f(3,lxInlet,hight1)=f(5,lxInlet,hight1);

f(6,lxInlet,hight1)=f(8,lxInlet,hight1);

f(7,lxInlet,hight1)=f(9,lxInlet,hight1);

%%%% Bottom Second Obstacle wall

f(3,lxOutlet,hight2)=f(5,lxOutlet,hight2);

f(6,lxOutlet,hight2)=f(8,lxOutlet,hight2);

f(7,lxOutlet,hight2)=f(9,lxOutlet,hight2);

%%%% Left Cavity wall

f(2,lenght1,lyLeftCavity)=f(4,lenght1,lyLeftCavity);

f(6,lenght1,lyLeftCavity)=f(8,lenght1,lyLeftCavity);

f(9,lenght1,lyLeftCavity)=f(7,lenght1,lyLeftCavity);

%%%% Right Cavity Wall

f(4,lenght1+lenght2,lyRightCavity)=f(2,lenght1+lenght2,lyRightCavity);

f(8,lenght1+lenght2,lyRightCavity)=f(6,lenght1+lenght2,lyRightCavity);

f(7,lenght1+lenght2,lyRightCavity)=f(9,lenght1+lenght2,lyRightCavity);

%%%% Bottom Cavity Wall

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

f(6,lxBottomCavity,1)=f(8,lxBottomCavity,1);

f(7,lxBottomCavity,1)=f(9,lxBottomCavity,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%

%%%%%%%%% rho and Velocity %%%%%%%%%%%%%%%%%%%%%%%%%%%

rho = reshape(sum(f),lx,ly);

u(:, = reshape(ex * reshape(f,9,lx*ly),lx,ly)./rho;

v(:, = reshape(ey * reshape(f,9,lx*ly),lx,ly)./rho;

u(1,lyInlet) = Uo;u(lxInlet,lyLeftCavity)= 0;u(lxOutlet,lyRightCavity)= 0;

v(1,lyInlet) = 0 ;v(lxInlet,lyLeftCavity)= 0;v(lxOutlet,lyRightCavity)= 0;

%

% v(lx,lyOutlet)=0;

% u(lx,lyOutlet)=u(lx-1,lyOutlet);

figure(2);

Uaverage=sqrt(u(:,:).^2+v(:,:).^2);

imagesc(Uaverage(:,ly:-1:1)’./Uo);

drawnow

int2str(t)

end