function [ux,uy,press,time,vel,velDifference,i] = anb(tmax,density,accel,omega,rrey,obst)
%
% ANB was originally written in FORTRAN
% LBGK Lattice Boltzmann, Bhatnagar-Gross-Krook
%
% Flow between two plates with obstacles
% The file par.m contains the parameters for the obstacles definitions
%
% Usage:
% [ux,uy,press,time,vel,velDifference,i] = anb(tmax,density,accel,omega,rrey,obst)
% where:
% ux = matrix of x-velocity components
% uy = matrix of y-velocity components
% press = matrix of pressure values
% i = number of time steps used
lx=400;
ly=200;
tmax=3000;
density=0.1;
accel=0.005;
omega=1/1.85;
rrey=10;
obst=zeros(lx,ly);
obst(:,[1,ly])=1;
% Node matrix (distribution matrix)
node = zeros(lx,ly,9);
nhlp = zeros(lx,ly,9);
vel = 0;
function [time, vel] = writeVelocity(lx,ly,time,obst,node,vel)
% obtains the macroscopic velocity
x = lx/2;
nfree = 0;
ux = 0;
for y = 1:ly
if( obst(x,y) == 0 )
dloc = 0;
for i = 1:9
dloc = dloc + node(x,y,i);
end
ux(x,y) = ( node(x,y,1) + node(x,y,5) + node(x,y,8)...
-( node(x,y,3) + node(x,y,6) + node(x,y,7) ) )/dloc;
nfree = nfree + 1;
end
end
vel = ux/nfree;
function nhlp = propagate(lx,ly,node,nhlp)
% nhlp = output matrix
% moves the fluid particles to the next lattice node
for x = 1:lx
for y = 1:ly
% Define north, south, east and west
yn = mod(y,ly) + 1;
xe = mod(x,lx) + 1;
ys = ly - mod(ly+1-y, ly);
xw = lx - mod(lx+1-x, lx);
% density propagation
nhlp(x ,y ,9) = node(x,y,9);
nhlp(xe ,y ,1) = node(x,y,1);
nhlp(x ,yn ,2) = node(x,y,2);
nhlp(xw ,y ,3) = node(x,y,3);
nhlp(x ,ys ,4) = node(x,y,4);
nhlp(xe ,yn ,5) = node(x,y,5);
nhlp(xw ,yn ,6) = node(x,y,6);
nhlp(xw ,ys ,7) = node(x,y,7);
nhlp(xe ,ys ,8) = node(x,y,8);
end
end
function node = redistribute(lx,ly,obst,node,accel,density)
% redistributes the density in the first lattice column
% drives the motion by modifying the distribution in the first lattice
% column
% weighting factors
t1 = density*accel/9;
t2 = density*accel/36;
for y = 1:ly
if (obst(1,y) == 0 & node(1,y,3)-t1 >= 0 & node(1,y,6)-t2 >= 0 & ...
node(1,y,7)-t2 >= 0 )
node(1,y,1) = node(1,y,1) + t1;
node(1,y,3) = node(1,y,3) - t1;
node(1,y,5) = node(1,y,5) + t2;
node(1,y,6) = node(1,y,6) - t2;
node(1,y,7) = node(1,y,7) - t2;
node(1,y,8) = node(1,y,8) + t2;
end
end
function node = relaxation(density, omega, lx,ly,node,nhlp,obst)
% implements the BGK approximation.
% Finds local equilibrium distributions and relaxes the distributions
% towards equilibrium
% Weighting factors in 9-bit LBGK
t0 = 4/9; t1 = 1/9; t2 = 1/36;
% c = speed of sound
% csqu = c squared
csqu = 1/3;
% Main relaxation loop
for x = 1:lx
for y = 1:ly
if (obst(x,y)==0)
% initialize variable dloc
dloc = 0;
for i = 1:9
dloc = dloc + nhlp(x,y,i);
end
% ux and uy velocity components
ux = (nhlp(x,y,1) + nhlp(x,y,5) + nhlp(x,y,8) -...
(nhlp(x,y,3) + nhlp(x,y,6) + nhlp(x,y,7)))/ dloc;
uy = (nhlp(x,y,2) + nhlp(x,y,5) + nhlp(x,y,6) -...
(nhlp(x,y,4) + nhlp(x,y,7) + nhlp(x,y,8)))/ dloc;
% usqu = magnitude of the velocity
usqu = ux^2 + uy^2;
% n-velocity components (n = lattice node connection vectors)
un(1) = ux;
un(2) = uy;
un(3) = -ux;
un(4) = -uy;
un(5) = ux + uy;
un(6) = -ux + uy;
un(7) = -ux - uy;
un(8) = ux - uy;
% Equilibrium densities
% index 0 in the fortran code is index 9 here
nequ(9) = t0*dloc*(1 - usqu/(2*csqu));
nequ(1) = t1*dloc*(1 + un(1)/csqu + un(1)^2/(2*csqu^2) - usqu/(2*csqu));
nequ(2) = t1*dloc*(1 + un(2)/csqu + un(2)^2/(2*csqu^2) - usqu/(2*csqu));
nequ(3) = t1*dloc*(1 + un(3)/csqu + un(3)^2/(2*csqu^2) - usqu/(2*csqu));
nequ(4) = t1*dloc*(1 + un(4)/csqu + un(4)^2/(2*csqu^2) - usqu/(2*csqu));
nequ(5) = t2*dloc*(1 + un(5)/csqu + un(5)^2/(2*csqu^2) - usqu/(2*csqu));
nequ(6) = t2*dloc*(1 + un(6)/csqu + un(6)^2/(2*csqu^2) - usqu/(2*csqu));
nequ(7) = t2*dloc*(1 + un(7)/csqu + un(7)^2/(2*csqu^2) - usqu/(2*csqu));
nequ(8) = t2*dloc*(1 + un(8)/csqu + un(8)^2/(2*csqu^2) - usqu/(2*csqu));
%
% relaxation towards the equilibrium
for i= 1:9
node(x,y) = nhlp(x,y,i) + omega*(nequ(i) - nhlp(x,y,i));
end
end
end
end
function [ux,uy,press] = writeResults(lx,ly,obst,node,density)
% obtains the macroscopic velocity components ux and uy and the pressure
% press
csqu = 1/3;
% in the original fortran code the matrices ux, uy, press are written to a
% file. Here the matrices ux, uy and press are passed to the function anb.m
ux = zeros(lx,ly);
uy = zeros(lx,ly);
press = zeros(lx,ly);
for y = 1:ly
for x = 1:lx
if (obst(x,y) == 1)
ux(x,y) = 0;
uy(x,y) = 0;
press(x,y) = density*csqu;
else
dloc = 0;
for i = 1:9
dloc = dloc + node(x,y,i);
end
% velocity components
ux(x,y) = ( node(x,y,1) + node(x,y,5) + node(x,y,8)...
-( node(x,y,3) + node(x,y,6) + node(x,y,7) ) )/dloc;
uy(x,y) = ( node(x,y,2) + node(x,y,5) + node(x,y,6)...
-( node(x,y,4) + node(x,y,7) + node(x,y,8) ) )/dloc;
press(x,y) = dloc*csqu;
end
end
end
function node = initdensity(lx,ly,density)
% initdensity(lx,ly,density,node) finds the initial density and a node-matrix
% the function is run before any time steps and sets the initial
% distributions
t0 = density*4/9;
t1 = density/9;
t2 = density/36;
for x = 1:lx
for y = 1:ly
node(x,y,9) = t0;
node(x,y,1) = t1;
node(x,y,2) = t1;
node(x,y,3) = t1;
node(x,y,4) = t1;
node(x,y,5) = t2;
node(x,y,6) = t2;
node(x,y,7) = t2;
node(x,y,8) = t2;
end
end
% Initial densities
node = initdensity(lx,ly,density);
%
% An iterative time step follows
% write interim results
[ux,uy,press] = writeResults(lx,ly,obst,node,density);
%
i = 1;
velDifference(1) = 1;
while i<tmax & velDifference(i) > 1/200000
% while i<tmax
i = i + 1;
if i >= 10 & mod(i, tmax/10) == 0
totalMass = checkDensity(lx,ly,node,i);
disp(sprintf('Time step: %d Total density: %0.5g', i, totalMass))
% disp(sprintf('Node difference: %0.5g', nodeDiff(node, nhlp)))
end
%disp('Step'); disp(i);
% A Control printout to find how long it takes to reach stationary flow
% Compare the old velocity matrices to the next velocity matrices
oldUx = ux; oldUy = uy;
% end of control printout
node = redistribute(lx,ly,obst,node,accel,density);
nhlp = propagate(lx,ly,node,nhlp);
for i=1:9
nhlp(i,obst)=node(opp(i),bbregion);
end
node = bounceback(lx,ly,obst,node,nhlp);
node = relaxation(density, omega, lx,ly,node,nhlp,obst);
% The difference in velocities between the two time steps
[ux,uy,press] = writeResults(lx,ly,obst,node,density);
velDifference(i) = sum(abs(ux - oldUx) + abs(uy - oldUy)) / sum(abs(ux) + abs(uy));
[time, vel] = writeVelocity(lx,ly,i,obst,node,vel);
end %while
[ux,uy,press] =
