Hi,
I have written a simple code for simulation of flow around a spherical object. I have followed guidelines of the example cylinder.m code given on the Palabos wiki site. Please let me know what I’m doing wrong:
%LB-3D Wind Tunnel
clear
clc
% y
% |
% |_____ x
% /
% /
% z
lx = 100; % number of cells in x-direction
ly = 50; % number of cells in y-direction
lz = 50; % number of cells in z-direction
%For sphere
obst_x = lx/5+1;
obst_y = ly/2;
obst_z = lz/2;
obst_r= ly/10+1;
%NOTE:no need for poiseulle flow, keep it for now
uMax = 0.1;
Re = 100; % Reynolds number
nu = uMax * 2.*obst_r / Re; % kinematic viscosity
omega = 1. / ( 3*nu+1./2.); % relaxation parameter
maxT = 400000; % total number of iterations
tPlot = 2; % cycles
%
%D3Q19 Lattice Constants
t= [1/3, 1/18, 1/36, 1/18, 1/36, 1/18, 1/36, 1/18, 1/36, 1/36, 1/18, 1/36, 1/36, 1/18, 1/36, 1/36, 1/36, 1/36, 1/36];
cx= [0, 1, 1, 0,-1,-1,-1, 0, 1, 0, 0, 0, 0, 0, 0, 1,-1,-1, 1];
cy= [0, 0, 1, 1, 1, 0,-1,-1,-1, 1, 0,-1,-1, 0, 1, 0, 0, 0, 0];
cz= [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1,-1,-1,-1, 1, 1,-1,-1];
opp=[1,6,7,8,9,2,3,4,5,13,14,15,10,11,12,18,19,16,17];
%
coly = [2:(ly-1)];
colz = [2:(lz-1)];
in=1;
out=lx;
[z,y,x] = meshgrid(1:ly,1:lx,1:lz);
%
%spherical obstacle for test, replace later with voxels
%
obst = ((x-obst_x).^2 + (y-obst_y).^2 + (z-obst_z).^2 <= obst_r.^2);
%obst(:,[1,ly],:)=1; % location of top and bottom
bbRegion = find(obst); %boolean mask for bounce back
% INITIAL CONDITION: Poiseuille profile at equilibrium
L = ly-2; y_phys = y-1.5;
ux = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys);
uy = zeros(lx,ly,lz);
uz = zeros(lx,ly,lz);
%
rho = 1;
for i=1:19
cu = 3*(cx(i)*ux + cy(i)*uy + cz(i)*uz);
fIn(i,:,:,:) = rho .* t(i) .*( 1 + cu + 1/2*(cu.*cu) - 3/2*(ux.^2 + uy.^2 + uz.^2) );
end
%
%%
% MAIN LOOP (TIME CYCLES)
for cycle = 1:maxT
% MACROSCOPIC VARIABLES
rho = sum(fIn);
ux = reshape ( (cx * reshape(fIn,19,lx*ly*lz)), 1,lx,ly,lz) ./rho;
uy = reshape ( (cy * reshape(fIn,19,lx*ly*lz)), 1,lx,ly,lz) ./rho;
uz = reshape ( (cz * reshape(fIn,19,lx*ly*lz)), 1,lx,ly,lz) ./rho;
% MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS
% Inlet: Poiseuille profile
y_phys = repmat((coly-1.5),lz-2,1);
%size(y_phys)
%size(ux(:,in,coly,colz))
ux(:,in,coly,colz) = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys);
uy(:,in,coly,colz) = 0;
rho(:,in,coly,colz) = 1 ./ (1-ux(:,in,coly,colz)) .* ( ...
sum(fIn([1,4,8,10,11,12,13,14,15],in,coly,colz)) + 2*sum(fIn([5,6,7,17,18],in,coly,colz)) );
% Outlet: Constant pressure
rho(:,out,coly,colz) = 1;
ux(:,out,coly,colz) = -1 + 1 ./ (rho(:,out,coly,colz)) .* ( ...
sum(fIn([1,4,8,10,11,12,13,14,15],out,coly,colz)) + 2*sum(fIn([9,13,16,19],out,coly,colz)) );
uy(:,out,coly,colz) = 0;
% MICROSCOPIC BOUNDARY CONDITIONS: INLET (Zou/He BC)
fIn(2,in,coly,colz) = fIn(6,in,coly,colz) + 1/3*rho(:,in,coly,colz).*ux(:,in,coly,colz);
fIn(3,in,coly,colz) = fIn(7,in,coly,colz) + 1/4 * ( fIn(8,in,coly,colz) - fIn(4,in,coly,colz) ) + 1/6*rho(:,in,coly,colz).*ux(:,in,coly,colz);
fIn(9,in,coly,colz) = fIn(5,in,coly,colz) + 1/4 * ( fIn(4,in,coly,colz) - fIn(8,in,coly,colz) ) + 1/6*rho(:,in,coly,colz).*ux(:,in,coly,colz);
fIn(16,in,coly,colz) = fIn(18,in,coly,colz) + 1/4 * ( fIn(14,in,coly,colz) - fIn(11,in,coly,colz) ) + 1/6*rho(:,in,coly,colz).*ux(:,in,coly,colz);
fIn(19,in,coly,colz) = fIn(17,in,coly,colz) + 1/4 * ( fIn(11,in,coly,colz) - fIn(14,in,coly,colz) ) + 1/6*rho(:,in,coly,colz).*ux(:,in,coly,colz);
% MICROSCOPIC BOUNDARY CONDITIONS: OULET (Zou/He BC)
fIn(6,out,coly,colz) = fIn(2,out,coly,colz) - 1/3*rho(:,out,coly,colz).*ux(:,out,coly,colz);
fIn(7,out,coly,colz) = fIn(3,out,coly,colz) - 1/4 * ( fIn(8,out,coly,colz) - fIn(4,out,coly,colz) ) - 1/6*rho(:,in,coly,colz).*ux(:,out,coly,colz);
fIn(5,out,coly,colz) = fIn(9,out,coly,colz) - 1/4 * ( fIn(4,out,coly,colz) - fIn(8,out,coly,colz) ) - 1/6*rho(:,in,coly,colz).*ux(:,out,coly,colz);
fIn(18,out,coly,colz) = fIn(16,out,coly,colz) - 1/4 * ( fIn(14,out,coly,colz) - fIn(11,out,coly,colz) ) - 1/6*rho(:,in,coly,colz).*ux(:,out,coly,colz);
fIn(17,out,coly,colz) = fIn(19,out,coly,colz) - 1/4 * ( fIn(11,out,coly,colz) - fIn(14,out,coly,colz) ) - 1/6*rho(:,in,coly,colz).*ux(:,out,coly,colz);
% COLLISION STEP
for i=1:19
cu = 3*(cx(i)*ux + cy(i)*uy + cz(i)*uz);
fOut(i,:,:,:) = rho .* t(i) .*( 1 + cu + 1/2*(cu.*cu) - 3/2*(ux.^2 + uy.^2 + uz.^2) );
% fOut(i,:,:,:) = fIn(i,:,:,:) - omega .* (fIn(i,:,:,:)-fEq(i,:,:,:));
end
% OBSTACLE (BOUNCE-BACK)
for i=1:19
fOut(i,bbRegion) = fIn(opp(i),bbRegion);
end
% STREAMING STEP
for i=1:19
fIn(i,:,:,:) = circshift(fOut(i,:,:,:), [0,cx(i),cy(i),cz(i)]);
end
% VISUALIZATION
if (mod(cycle,tPlot)==1)
u = reshape(sqrt(ux.^2+uy.^2 + uz.^2),lx,ly,lz);
u(bbRegion) = nan;
imagesc((u(:,:,lz/2))');
axis equal off; drawnow
end
end