Hi Everyone,
I’m a new to the world of LB simulations with background in multiphysics modelling (namely COMSOL). I’m trying to simulate a simple droplet spreading on a flat surface using a D2Q9 Lattice and a BGK model in MATLAB. The m-files by Adriano Sciacovelli and Orestis Malaspinas et al. have been great help in starting me up but I’m finding some difficulties in implementing the adhesion force as described in Sukop’s Introduction to LBM book.
I originally suspected that the mistake in the code came from using array manipulations with the adhesion ‘switch’ implementation instead of using for-loop statements but then I didn’t have much success when coding the latter. In fact for some reason the density wasn’t conserved in the system and continually increased!
Anyways here is the code for a SCMP simulation (where the mass is at least conserved); it works fine if s(20:80,52:55) = 1 is commented out. There are 2 initial density conditions (rho0) to play with. I hope someone could give me a hint about how to fix it.
Thank you very much in advance 
Laith
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
clf
% Model constants
Nx = 101;
Ny = 101;
rho0=200; %liquid density [kg/m^3]
tau = 1; %BGK relaxation coeff [ts].
nu =(1/3)(tau-0.5);%kinematic viscosity in lattice units [lu^2/ts]
omega=1/tau;
lb = Ny/5; %LBM length scale of lbp [lu];
Nt = 1000; %time steps in [ts]
RT = 1/3; %value of RT in D2Q9;
G = -120; %interaction constant for SCMP
Gad=-150; %adsorption coefficient for SCMP
psi0 = 4;
%D2Q9 lattice constants
w = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
ex = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
ey = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
opp = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];
%descritize domains
[x,y] = meshgrid(1:Nx,1:Ny);
% initialize domain
s = zeros(Nx,Ny);
%s(20:80,51:55) = 1;
obst=find(s==1);
ux = zeros(Nx,Ny);
uy = ux;
rho = rho0+0.1*(rand(Nx,Ny)>0.55);
%rho = 81+(540-81)(lb.^2>((y-50).^2+1(x-30).^2));
fIn = zeros(9,Nx,Ny);
fOut = zeros(9,Nx,Ny);
fEq = zeros(9,Nx,Ny);
uxtemp = zeros(9,Nx,Ny);
uytemp = zeros(9,Nx,Ny);
Fxtemp = zeros(9,Nx,Ny);
Fytemp = Fxtemp;
Fadxtemp = zeros(9,Nx,Ny);
Fadytemp = Fadxtemp;
for i=1:9
eu = ex(i).*ux + ey(i).*uy;
fIn(i,:, = w(i).rho.(1);%+3.eu+(9/2)(eu.eu)-(3/2)(ux.^2+uy.^2));
end
% main loop
for count=1:Nt
% calculate Macroscopic properties
rho = sum(fIn);
psi = psi0.*exp(-rho0./rho);
psi(1,obst) = 0; %I'm not sure if this is needed
P = (rho./3) + (G/6)*psi.^2;
%calculate forcing based on long range attractive (cohesive) forces
for i=1:9
if i>1
Fxtemp(i,:,:) = ex(i).*circshift(w(i).*psi(1,:,:),[0,ex(i),ey(i)]);
Fytemp(i,:,:) = ey(i).*circshift(w(i).*psi(1,:,:),[0,ex(i),ey(i)]);
end
Fadxtemp(i,:,:) = ex(i).*circshift(w(i).*s,[ex(i),ey(i)]);
Fadytemp(i,:,:) = ey(i).*circshift(w(i).*s,[ex(i),ey(i)]);
uxtemp(i,:,:) = fIn(i,:,:).*ex(i);
uytemp(i,:,:) = fIn(i,:,:).*ey(i);
end
%calculate force
Fx = -G*psi.*sum(Fxtemp);
Fy = -G*psi.*sum(Fytemp);
Fadx = -Gad*psi.*sum(Fadxtemp);
Fady = -Gad*psi.*sum(Fadytemp);
%calculate velocity
ux = -(Fadx.*tau./rho)-(Fx.*tau./rho) + sum(uxtemp)./rho;
uy = -(Fady.*tau./rho)-(Fy.*tau./rho) + sum(uytemp)./rho;
% apply periodic BC
%(circshift should give a periodic BC?)
% collide
for i=1:9
eu = ex(i).*ux + ey(i).*uy;
fEq(i,:,:) = w(i).*rho.*(1+3.*eu+(9/2)*(eu.*eu)-(3/2)*(ux.^2+uy.^2));
fOut(i,:,:) = fIn(i,:,:) - omega.*(fIn(i,:,:)-fEq(i,:,:));
end
% stream
for i=1:9
fIn(i,:,:) = circshift(fOut(i,:,:), [0,ex(i),ey(i)]);
end
% plots
u = reshape(sqrt(ux.^2+uy.^2),Nx,Ny);
rhop = reshape(rho,Nx,Ny);
imagesc(rhop'); axis equal off; drawnow;
rsum(count) = sum(sum(rhop));
%plot(1:count,rsum);drawnow
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%