Simple LBM code doubt (code included)

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