Dear All,
Problem solved - I have had a typo in the source. I thank to all that have had a look at the post.
I have taken as a template the matlab script (shanchen.m) of Orestis Malaspinas, Andrea Parmigiani and Jonas Latt that is posted on this site and modified it to run in 3D using the D3Q19 lattice according to “Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models”, Huang, Thorne, Schaap, Sukop, PHYSICAL REVIEW E 76, 066701, 2007
The modified script seems to run OK, for small amplitudes of the mollecular interaction force G (say G= -0.5), but diverges for larger ones (say G= -0.9). For small G a diffuse interface is obtained and larger Gs do not produce a stable solution.
In contrast the 2D version shanchen.m is much more stable and sharp interfaces and stable solutions could be obtained for appropriate values of G.
I doubt that such a difference in behaviour of the 2D and 3D implementations have any theoretical origin, but any comments on the issue are greatly appreciated.
I have copied the modified script below.
Kind regards,
Yulian
function [r1, r2]= shanchen_3D(maxT)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% shanChen.m: Multi-component fluid, using a LB method,
% based on the Shan-Chen model
% [X.Shan and H.Chen, http://dx.doi.org/10.1103/PhysRevE.47.1815].
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lattice Boltzmann sample, written in Matlab
% Copyright © 2008 Orestis Malaspinas, Andrea Parmigiani, Jonas Latt
% Address: EPFL-STI-LIN Station 9
% E-mail: orestis.malaspinas@epfl.ch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
% You should have received a copy of the GNU General Public
% License along with this program; if not, write to the Free
% Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
% Boston, MA 02110-1301, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%clear
%clf
% GENERAL FLOW CONSTANTS
ly = 21;
lx = 21;
lz = 21;
G = -0.9; % Amplitude of the molecular interaction force
nu1= 1/6;
nu2= 1/6;
omega1 = 1/(3nu1+0.5); % Relaxation parameter for fluid 1
omega2 = 1/(3nu2+0.5); % Relaxation parameter for fluid 2
%maxT = 3000; % total number of iterations
tPlot = 1; % iterations between successive graphical outputs
out_count= 0;
r1= zeros(round(maxT/tPlot)+1, lx, ly, lz);
r2= zeros(round(maxT/tPlot)+1, lx, ly, lz);
%% D2Q9 LATTICE CONSTANTS
%tNS = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
%cxNS = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
%cyNS = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
%oppNS = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];
% D3Q19 LATTICE CONSTANTS
tNS = [1/3, 1/18,1/18,1/18,1/18,1/18,1/18, 1/36,1/36,1/36,1/36,1/36,1/36,1/36,1/36,1/36,1/36,1/36,1/36];
cxNS = [ 0, 1,-1, 0, 0, 0, 0, 1, 1,-1,-1, 1,-1, 1,-1, 0, 0, 0, 0];
cyNS = [ 0, 0, 0, 1,-1, 0, 0, 1,-1, 1,-1, 0 ,0, 0, 0, 1, 1,-1,-1];
czNS = [ 0, 0, 0, 0, 0, 1,-1, 0, 0, 0, 0, 1, 1,-1,-1, 1,-1 ,1,-1];
[y,x,z] = meshgrid(1:ly,1:lx,lz);
drho = 0.001;
delta_rho = -drho*(1-2.0*rand(lx,ly,lz));
cent_x1= round(lx/3); cent_y1= round(ly/2); cent_z1= round(lz/2);
cent_x2= round(2*lx/3); cent_y2= round(ly/2); cent_z2= round(lz/2);
rad= round(lx/10);
circ_in= zeros(lx,ly,lz); circ_out= ones(lx,ly,lz);
circ_in(cent_x1-rad:cent_x1+rad,cent_y1-rad:cent_y1+rad,cent_z1-rad:cent_z1+rad)= 1;
circ_in(cent_x2-rad:cent_x2+rad,cent_y2-rad:cent_y2+rad,cent_z2-rad:cent_z2+rad)= 1;
circ_out(cent_x1-rad:cent_x1+rad,cent_y1-rad:cent_y1+rad,cent_z1-rad:cent_z1+rad)= 0;
circ_out(cent_x2-rad:cent_x2+rad,cent_y2-rad:cent_y2+rad,cent_z2-rad:cent_z2+rad)= 0;
% INITIAL CONDITION FOR BOTH DISTRIBUTION FUNCTIONS: (T=0) ==> TIn(i) = t(i)
for i=1:19
% fIn(i,1:lx,1:ly,1:lz) = tNS(i)circ_in2;
% gIn(i,1:lx,1:ly,1:lz) = tNS(i)circ_out2;
fIn(i,1:lx,1:ly,1:lz) = tNS(i).(1.0 + delta_rho);
gIn(i,1:lx,1:ly,1:lz) = tNS(i).(1.0 - delta_rho);
end
rho1 = reshape(sum(fIn),lx,ly,lz);
% MAIN LOOP (TIME CYCLES)
Gomega1 = G/omega1;
Gomega2 = G/omega2;
for cycle = 1:maxT
cycle
% MACROSCOPIC VARIABLES
rho1 = sum(fIn);
rho2 = sum(gIn);
jx1 = reshape ( (cxNS * reshape(fIn,19,lxlylz)), 1,lx,ly,lz);
jy1 = reshape ( (cyNS * reshape(fIn,19,lxlylz)), 1,lx,ly,lz);
jz1 = reshape ( (czNS * reshape(fIn,19,lxlylz)), 1,lx,ly,lz);
jx2 = reshape ( (cxNS * reshape(gIn,19,lx*ly*lz)), 1,lx,ly,lz);
jy2 = reshape ( (cyNS * reshape(gIn,19,lx*ly*lz)), 1,lx,ly,lz);
jz2 = reshape ( (czNS * reshape(gIn,19,lx*ly*lz)), 1,lx,ly,lz);
rhoTot_OMEGA = rho1*omega1 + rho2*omega2;
uTotX = (jx1*omega1+jx2*omega2) ./ rhoTot_OMEGA;
uTotY = (jy1*omega1+jy2*omega2) ./ rhoTot_OMEGA;
uTotZ = (jz1*omega1+jz2*omega2) ./ rhoTot_OMEGA;
rhoContrib1x = 0.0;
rhoContrib2x = 0.0;
rhoContrib1y = 0.0;
rhoContrib2y = 0.0;
rhoContrib1z = 0.0;
rhoContrib2z = 0.0;
for i=2:19
temp1CS= circshift(rho1*tNS(i), [0,cxNS(i),cyNS(i),czNS(i)]);
rhoContrib1x = rhoContrib1x + temp1CS*cxNS(i);
rhoContrib1y = rhoContrib1y + temp1CS*cyNS(i);
rhoContrib1z = rhoContrib1z + temp1CS*czNS(i);
% rhoContrib1x = rhoContrib1x + circshift(rho1*tNS(i), [0,cxNS(i),cyNS(i),czNS(i)])cxNS(i);
% rhoContrib1y = rhoContrib1y + circshift(rho1tNS(i), [0,cxNS(i),cyNS(i),czNS(i)])cyNS(i);
% rhoContrib1z = rhoContrib1z + circshift(rho1tNS(i), [0,cxNS(i),cyNS(i),czNS(i)])*czNS(i);
temp1CS= circshift(rho2*tNS(i), [0,cxNS(i),cyNS(i),czNS(i)]);
rhoContrib2x = rhoContrib2x + temp1CS*cxNS(i);
rhoContrib2y = rhoContrib2y + temp1CS*cyNS(i);
rhoContrib2z = rhoContrib2z + temp1CS*czNS(i);
% rhoContrib2x = rhoContrib2x + circshift(rho2*tNS(i), [0,cxNS(i),cyNS(i),czNS(i)])cxNS(i);
% rhoContrib2y = rhoContrib2y + circshift(rho2tNS(i), [0,cxNS(i),cyNS(i),czNS(i)])cyNS(i);
% rhoContrib2z = rhoContrib2z + circshift(rho2tNS(i), [0,cxNS(i),cyNS(i),czNS(i)])*czNS(i);
end
uTotX1 = uTotX - Gomega1.*rhoContrib2x; %POTENTIAL CONTRIBUTION OF FLUID 2 ON 1
uTotY1 = uTotY - Gomega1.*rhoContrib2y;
uTotZ1 = uTotZ - Gomega1.*rhoContrib2z;
uTotX2 = uTotX - Gomega2.*rhoContrib1x; %POTENTIAL CONTRIBUTION OF FLUID 2 ON 1
uTotY2 = uTotY - Gomega2.*rhoContrib1y;
uTotZ2 = uTotZ - Gomega2.*rhoContrib1z;
% COLLISION STEP FLUID 1 AND 2
for i=1:19
cuNS1 = 3*(cxNS(i)*uTotX1+cyNS(i)*uTotY1+czNS(i)uTotZ1);
cuNS2 = 3(cxNS(i)*uTotX2+cyNS(i)*uTotY2+czNS(i)*uTotZ2);
fEq(i,:,:,:) = rho1 .* tNS(i) .* ...
( 1 + cuNS1 + 0.5*(cuNS1.*cuNS1) - 1.5*(uTotX1.^2+uTotY1.^2) );
gEq(i,:,:,:) = rho2 .* tNS(i) .* ...
( 1 + cuNS2 + 0.5*(cuNS2.*cuNS2) - 1.5*(uTotX2.^2+uTotY2.^2) );
fOut(i,:,:,:) = fIn(i,:,:,:) - omega1 .* (fIn(i,:,:,:)-fEq(i,:,:,:));
gOut(i,:,:,:) = gIn(i,:,:,:) - omega2 .* (gIn(i,:,:,:)-gEq(i,:,:,:));
end
% STREAMING STEP FLUID 1 AND 2
for i=1:19
fIn(i,:,:, = circshift(fOut(i,:,:,:), [0,cxNS(i),cyNS(i),czNS(i)]);
gIn(i,:,:, = circshift(gOut(i,:,:,:), [0,cxNS(i),cyNS(i),czNS(i)]);
end
% VISUALIZATION
if(mod(cycle,tPlot)==0)
out_count= out_count+ 1;
rho1 = reshape(rho1,lx,ly,lz);
rho2 = reshape(rho2,lx,ly,lz);
r1(out_count,:,:,:)= rho1;
r2(out_count,:,:,:)= rho2;
subplot(2,2,1);
imagesc(rho1(:,:,round(lz/2)')); colorbar
sum_rho1=sum(rho1(:));
title(['rho1= ', num2str(sum_rho1)])
axis equal off; drawnow
subplot(2,2,2);
imagesc(rho2(:,:,round(lz/2)')); colorbar
sum_rho2=sum(rho2(:));
title(['rho2= ', num2str(sum_rho2)])
axis equal off; drawnow
subplot(2,2,3);
line(cycle,sum_rho1,'marker','.', 'markeredgecolor','r');
subplot(2,2,4);
line(cycle,sum_rho2,'marker','*','markeredgecolor','b');
pause(0.2)
end
end