# Simulating Contact angle on solids with LBM

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 R
T in D2Q9;
G = -120; %interaction constant 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;

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

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);

%calculate velocity

% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

this code:

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

be changed by
eu=3*(ex(i).*ux + ey(i).*uy);

is it right?