%3D Lattice Boltzmann (BGK) model of a fluid.
% D3Q19 model. At each timestep, particle densities propagate
% outwards in the directions indicated in the figure. An
% equivalent 'equilibrium' density is found, and the densities
% relax towards that state, in a proportion governed by omega.
% Iain Haslam, March 2006.
nx=12;ny=nx;nz=nx; omega=1.0; density=1.0;t1=1/3; t2=1/18; t3=1/36;
F=repmat(density/19,[nx ny nz 19]); FEQ=F; matsize=nx*ny*nz;
CI=[0:matsize:matsize*19];
BOUND=zeros(nx,ny,nz);
for i=1:nx, for j=1:ny, for k=1:nz
BOUND(i,j,k)=((i-5)^2+(j-6)^2+(k-7)^2)<6;
end, end, end
BOUND(:,:,1)=1;BOUND(:,1,:)=1;
ON=find(BOUND); %matrix offset of each Occupied Node
TO_REFLECT=[ON+CI(2) ON+CI(3) ON+CI(4) ON+CI(5) ON+CI(6) ON+CI(7) ON+CI(8) ...
ON+CI(9) ON+CI(10) ON+CI(11) ON+CI(12) ON+CI(13) ON+CI(14) ON+CI(15) ...
ON+CI(16) ON+CI(17) ON+CI(18) ON+CI(19)];
REFLECTED=[ON+CI(3) ON+CI(2) ON+CI(5) ON+CI(4) ON+CI(7) ON+CI(6) ON+CI(11) ...
ON+CI(10) ON+CI(9) ON+CI(8) ON+CI(15) ON+CI(14) ON+CI(13) ON+CI(12) ...
ON+CI(19) ON+CI(18) ON+CI(17) ON+CI(16)];
avu=1; prevavu=1; ts=0; deltaU=1e-7; numactivenodes=sum(sum(sum(1-BOUND)));
while (ts<4000 & 1e-10<abs((prevavu-avu)/avu)) | ts<100
% Propagate
%nearest-neighbours
F(:,:,:,2)=F(:,:,[nz 1:nz-1],2);
F(:,:,:,3)=F(:,:,[2:nz 1],3);
F(:,:,:,4)=F(:,[ny 1:ny-1],:,4);
F(:,:,:,5)=F(:,[2:ny 1],:,5);
F(:,:,:,6)=F([nx 1:nx-1],:,:,6);
F(:,:,:,7)=F([2:nx 1],:,:,7);
%next-nearest neighbours
F(:,:,:,8)= F([nx 1:nx-1],[ny 1:ny-1],:,8);
F(:,:,:,9)= F([nx 1:nx-1],[2:ny 1],:,9);
F(:,:,:,10)=F([2:nx 1],[ny 1:ny-1],:,10);
F(:,:,:,11)=F([2:nx 1],[2:ny 1],:,11);
F(:,:,:,12)=F([nx 1:nx-1],:,[nz 1:nz-1],12);
F(:,:,:,13)=F([nx 1:nx-1],:,[2:nz 1],13);
F(:,:,:,14)=F([2:nx 1],:,[nz 1:nz-1],14);
F(:,:,:,15)=F([2:nx 1],:,[2:nz 1],15);
F(:,:,:,16)=F(:,[ny 1:ny-1],[nz 1:nz-1],16);
F(:,:,:,17)=F(:,[ny 1:ny-1],[2:nz 1],17);
F(:,:,:,18)=F(:,[2:ny 1],[nz 1:nz-1],18);
F(:,:,:,19)=F(:,[2:ny 1],[2:nz 1],19);
评论0