Solved! Shan Chen in 3D (following the matlab source by Malaspinas, Parmigiani and Latt)

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/(3
nu2+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(rho1
tNS(i), [0,cxNS(i),cyNS(i),czNS(i)])cyNS(i);
% rhoContrib1z = rhoContrib1z + circshift(rho1
tNS(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(rho2
tNS(i), [0,cxNS(i),cyNS(i),czNS(i)])cyNS(i);
% rhoContrib2z = rhoContrib2z + circshift(rho2
tNS(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,:,:,:slight_smile: = circshift(fOut(i,:,:,:), [0,cxNS(i),cyNS(i),czNS(i)]);
gIn(i,:,:,:slight_smile: = 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

Dear Yulian ,

I am a new person begin to study the LBM, and not very familiar with the Matlab.

I was wondering the how can I get a figure (like Fig.3 in the Huang’s paper you mentioned) through Matlab?

Thank you for your help!

Best,
Iris