MATLAB code: Lattice Boltzmann method (D2Q9) for fluid flow passing through obstacles
Date:
%Main code: Main.m
clear all;
nx = 501; ny = 81;
f = zeros(nx,ny,9); feq = zeros(nx,ny,9);
u = zeros(nx,ny); v = zeros(nx,ny);
rho = ones(nx,ny); x = zeros(nx); y = zeros(ny);
Tm = zeros(nx); w(9) = zeros; Tvm = zeros(nx);
w = [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];
c2 = 1.0/3.0;
dx = 1.0; dy = 1.0;
xl = (nx-1)/(ny-1); yl = 1.0;
dx = xl/(nx-1);
dy = yl/(ny-1);
x = (0:dx:xl);
y = (0:dy:yl);
uo = 0.1;
alpha = 0.04;
Re = uo*(ny-1)/alpha;
omega = 1.0/(3.0*alpha+0.5);
count = 0; tol = 1.0e-4; error = 10.0; erso = 0.0;
%setting velocity
for j = 2:ny-1
u(1,j) = uo;
end
%Main Loop
while error > tol
%Collitions
[f] = collition(nx,ny,u,v,cx,cy,omega,f,rho,w);
%Streaming
[f] = stream(f);
%Boundary condition
[f] = boundary(nx,ny,f,uo,rho);
%Solid Obsticle
[f] = obstc(nx,ny,f,uo,rho);
%Calculate rho, u, v
[rho,u,v] = ruv(nx,ny,f);
count = count + 1;
ers = 0.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;
end
%Show result
result(nx,ny,x,y,u,v,uo,rho);
%Collition.m
function [f] = collition(nx,ny,u,v,cx,cy,omega,f,rho,w)
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)*w(k)*(1.0+3.0*t2+4.5*t2*t2-1.5*t1);
f(i,j,k) = (1.0-omega)*f(i,j,k) + omega*feq(i,j,k);
end
end
end
end
%boundary.m
function [f] = boundary(nx,ny,f,uo,rho)
nyl = (ny-1)/2;
%right boundary, outflow
for j = 1:ny
f(nx,j,3) = f(nx-1,j,3);
f(nx,j,7) = f(nx-1,j,7);
f(nx,j,6) = f(nx-1,j,6);
end
%bottom and top boundaries, bounce back
for i = 1:nx
f(i,1,2) = f(i,1,4);
f(i,1,5) = f(i,1,7);
f(i,1,6) = f(i,1,8);
f(i,ny,4) = f(i,ny,2);
f(i,ny,7) = f(i,ny,5);
f(i,ny,8) = f(i,ny,6);
u(i,1) = 0.0; v(i,1) = 0.0;
u(i,ny) = 0.0; v(i,ny) = 0.0;
end
%Left boundary, inflow
for j = 2:ny-1
f(1,j,1) = f(1,j,3) + 2.0*rho(1,j)*uo/3.0;
f(1,j,5) = f(1,j,7) - 0.5*(f(1,j,2)-f(1,j,4))+rho(1,j)*uo/6.0;
f(1,j,8) = f(1,j,6) + 0.5*(f(1,j,2)-f(1,j,4))+rho(1,j)*uo/6.0;
u(1,j) = uo; v(1,j) = 0.0;
end
end
%obstc.m
function [f] = obstc(nx,ny,f,uo,rho)
nxL1 = 50;
nxR1 = 70;
nyB = 20;
nyT = 60;
nxL2 = 120;
nxR2 = 140;
nxL3 = 190;
nxR3 = 210;
%For the boundary of obsticle, use bounce back
%for Solid 1
for i = nxL1:nxR1
f(i,nyT,2) = f(i,nyT,4);
f(i,nyT,5) = f(i,nyT,7);
f(i,nyT,6) = f(i,nyT,8);
end
for i = nxL1:nxR1
f(i,nyB,4) = f(i,nyT,2);
f(i,nyB,7) = f(i,nyT,5);
f(i,nyB,8) = f(i,nyT,6);
end
for j = nyB:nyT
f(nxR1,j,1) = f(nxR1,j,3);
f(nxR1,j,5) = f(nxR1,j,7);
f(nxR1,j,8) = f(nxR1,j,6);
end
for j = nyB:nyT
f(nxL1,j,3) = f(nxL1,j,1);
f(nxL1,j,7) = f(nxL1,j,5);
f(nxL1,j,6) = f(nxL1,j,8);
end
%for Solid 2
for i = nxL2:nxR2
f(i,nyT,2) = f(i,nyT,4);
f(i,nyT,5) = f(i,nyT,7);
f(i,nyT,6) = f(i,nyT,8);
end
for i = nxL2:nxR2
f(i,nyB,4) = f(i,nyT,2);
f(i,nyB,7) = f(i,nyT,5);
f(i,nyB,8) = f(i,nyT,6);
end
for j = nyB:nyT
f(nxR2,j,1) = f(nxR2,j,3);
f(nxR2,j,5) = f(nxR2,j,7);
f(nxR2,j,8) = f(nxR2,j,6);
end
for j = nyB:nyT
f(nxL2,j,3) = f(nxL2,j,1);
f(nxL2,j,7) = f(nxL2,j,5);
f(nxL2,j,6) = f(nxL2,j,8);
end
%for Solid 3
for i = nxL3:nxR3
f(i,nyT,2) = f(i,nyT,4);
f(i,nyT,5) = f(i,nyT,7);
f(i,nyT,6) = f(i,nyT,8);
end
for i = nxL3:nxR3
f(i,nyB,4) = f(i,nyT,2);
f(i,nyB,7) = f(i,nyT,5);
f(i,nyB,8) = f(i,nyT,6);
end
for j = nyB:nyT
f(nxR3,j,1) = f(nxR3,j,3);
f(nxR3,j,5) = f(nxR3,j,7);
f(nxR3,j,8) = f(nxR3,j,6);
end
for j = nyB:nyT
f(nxL3,j,3) = f(nxL3,j,1);
f(nxL3,j,7) = f(nxL3,j,5);
f(nxL3,j,6) = f(nxL3,j,8);
end
end
%result.m
function result(nx,ny,x,y,u,v,uo,rho)
figure(13)
nxL1 = 50;
nxR1 = 70;
nyB = 20;
nyT = 60;
nxL2 = 120;
nxR2 = 140;
nxL3 = 190;
nxR3 = 210;
for i = nxL1:nxR1
for j = nyB:nyT
u(i,j) = 0.0;
v(i,j) = 0.0;
end
end
for i = nxL2:nxR2
for j = nyB:nyT
u(i,j) = 0.0;
v(i,j) = 0.0;
end
end
for i = nxL3:nxR3
for j = nyB:nyT
u(i,j) = 0.0;
v(i,j) = 0.0;
end
end
subplot(2,1,1)
ll = linspace(0,500/80,501);ll2 = linspace(0,1,81);kk = 5; kk2 = 5;
[xxl,yyl] = meshgrid(ll,ll2); xxl = xxl’; yyl = yyl’;
quiver(xxl(1:kk:end,1:kk2:end),yyl(1:kk:end,1:kk2:end),1u(1:kk:end,1:kk2:end),1v(1:kk:end,1:kk2:end),0)
axis([0 500/80 0 1])
box on; axis image
patch([50/80 70/80 70/80 50/80],[20/80 20/80 60/80 60/80],’k’)
patch([120/80 120/80 140/80 140/80],[20/80 20/80 60/80 60/80],’k’)
patch([120/80 140/80 140/80 120/80],[20/80 20/80 60/80 60/80],’k’)
patch([190/80 210/80 210/80 190/80],[20/80 20/80 60/80 60/80],’k’)
subplot(2,1,2)
ll = linspace(0,500/80,501);ll2 = linspace(0,1,81);kk = 2; kk2 = 2;
[xxl,yyl] = meshgrid(ll,ll2); u = u’; v = v’;
streamslice(xxl(1:kk:end,1:kk2:end),yyl(1:kk:end,1:kk2:end),1u(1:kk:end,1:kk2:end),1v(1:kk:end,1:kk2:end))
axis([0 500/80 0 1])
box on; axis image
patch([50/80 70/80 70/80 50/80],[20/80 20/80 60/80 60/80],’k’)
patch([120/80 120/80 140/80 140/80],[20/80 20/80 60/80 60/80],’k’)
patch([120/80 140/80 140/80 120/80],[20/80 20/80 60/80 60/80],’k’)
patch([190/80 210/80 210/80 190/80],[20/80 20/80 60/80 60/80],’k’)
end
%ruv.m
function [rho,u,v] = ruv(nx,ny,f)
rho = sum(f,3);
%Calculate velocities
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
%stream.m
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