MATLAB code: Lattice Boltzmann method (D2Q9) for binary phase-field surfactant model
Date:
%Main code
clear all;
nx = 121; ny = 121;
f = zeros(nx,ny,9);
f2 = zeros(nx,ny,9);
feq = zeros(nx,ny,9);
u = zeros(nx+2,ny+2); %phase-field
v = zeros(nx+2,ny+2); %surfactant
mu = zeros(nx,ny); %chemical potential for phase-field
mv = zeros(nx,ny); %chemical potential for surfactant
x = zeros(nx);
y = zeros(ny);
w(9) = zeros;
w = [1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36 4/9];
cx = [1 0 -1 0 1 -1 -1 1 0];
cy = [0 1 0 -1 1 1 -1 -1 0];
c2 = 1.0/3.0;
xl = 120.0; yl = 120.0;
dx = xl/(nx-1);
dy = yl/(ny-1);
dt = 1.0;
x = (0:dx:xl);
y = (0:dy:yl);
T = 0.001;
Pi = 0.0001*1.35;
Ex = 0.1*0.11;
Pi = 1.35; %1.35;
Ex = 0.5; %0.5;
beta = sqrt(0.09/(8*4.5));
kk = 4.5*beta;
kb = 0.5betaPi;
W = 0.5*beta/Ex;
eta = 1.0;
eta2 = 1.0;
M = 0.5;
count = 0;
tau = M/(dtetac2)+0.5;
tau2 = M/(dteta2c2)+0.5;
omega = 1./tau;
omega2 = 1./tau2;
%setting initial condition
for i = 2:nx+1
for j = 2:ny+1
u(i,j) = 0.55 + 0.01*(2*rand-1);
v(i,j) = 0.2;
end
end
%Main Loop
while count < 3000
%Collitions
[f] = collition(nx,ny,u,mu,cx,cy,eta,omega,f,w);
[f2] = collition(nx,ny,v,mv,cx,cy,eta2,omega2,f2,w);
%Streaming
[f] = stream(f);
[f2] = stream(f2);
%Boundary condition
[f] = boundary(nx,ny,f);
[f2] = boundary(nx,ny,f2);
%Calculate u,mu
[u,mu,v,mv] = ruv(nx,ny,u,mu,v,mv,f,f2,w,kk,beta,kb,W,T,c2,dt);
count = count + 1
end
%Show result
result(nx,ny,x,y,u,v);
%boundary.m
function [f] = boundary(nx,ny,f)
%right boundary
for j = 1:ny
f(nx,j,3) = f(1,j,3);
f(nx,j,7) = f(1,j,7);
f(nx,j,6) = f(1,j,6);
end
%bottom and top boundaries
for i = 1:nx
f(i,1,2) = f(i,ny,2);
f(i,1,5) = f(i,ny,5);
f(i,1,6) = f(i,ny,6);
f(i,ny,4) = f(i,1,4);
f(i,ny,7) = f(i,1,7);
f(i,ny,8) = f(i,1,8);
end
%Left boundary
for j = 1:ny
f(1,j,1) = f(nx,j,1);
f(1,j,5) = f(nx,j,5);
f(1,j,8) = f(nx,j,8);
end
end
%collition.m
function [f] = collition(nx,ny,u,mu,cx,cy,eta,omega,f,w)
for j = 1:ny
for i = 1:nx
for k = 1:9
if(k == 9)
feq(i,j,k) = u(i+1,j+1) + (w(k)-1.0)*eta*mu(i,j);
else
feq(i,j,k) = w(k)*eta*mu(i,j);
end
f(i,j,k) = (1.0-omega)*f(i,j,k) + omega*feq(i,j,k);
end
end
end
end
%result.m
function result(nx,ny,x,y,u,v)
figure(15)
xl = linspace(0,100,nx);yl = xl;
[xxl,yyl]=meshgrid(xl,yl);
A = u’; B = v’;
subplot(1,2,1)
surf(xxl,yyl,A(2:nx+1,2:ny+1));
shading interp;
%colormap jet
axis([0 100 0 100])
box on; axis image
subplot(1,2,2)
surf(xxl,yyl,B(2:nx+1,2:ny+1));
shading interp;
%colormap jet
axis([0 100 0 100])
box on; axis image
end
%ruv.m
function [u,mu,v,mv] = ruv(nx,ny,u,mu,v,mv,f,f2,w,kk,beta,kb,W,T,c2,dt)
toll = 0.01; %1.0e-10;
etaa = 0.001;
Lapu = zeros(nx,ny);
Gv_x = zeros(nx,ny);
Gv_y = zeros(nx,ny);
Gu_x = zeros(nx,ny);
Gu_y = zeros(nx,ny);
u(2:nx+1,2:ny+1) = sum(f,3);
v(2:nx+1,2:ny+1) = sum(f2,3);
%periodic condition
u(1,:) = u(nx+1,:);
u(nx+2,:) = u(2,:);
u(:,1) = u(:,ny+1);
u(:,ny+2) = u(:,2);
%periodic condition
v(1,:) = v(nx+1,:);
v(nx+2,:) = v(2,:);
v(:,1) = v(:,ny+1);
v(:,ny+2) = v(:,2);
for i = 1:nx
for j = 1:ny
lapp = 0.0;
lapp = lapp + w(1)*( u(i+2,j+1) - 2*u(i+1,j+1) + u(i,j+1) );
lapp = lapp + w(5)*( u(i+2,j+2) - 2.0*u(i+1,j+1) +u(i,j) );
lapp = lapp + w(2)*( u(i+1,j+2) - 2.0*u(i+1,j+1) +u(i+1,j) );
lapp = lapp + w(6)*( u(i,j+2) - 2.0*u(i+1,j+1) + u(i+2,j) );
lapp = lapp + w(3)*( u(i,j+1) - 2.0*u(i+1,j+1) + u(i+2,j+1) );
lapp = lapp + w(7)*( u(i,j) - 2.0*u(i+1,j+1) + u(i+2,j+2) );
lapp = lapp + w(4)*( u(i+1,j) - 2.0*u(i+1,j+1) + u(i+1,j+2) );
lapp = lapp + w(8)*( u(i+2,j) - 2.0*u(i+1,j+1) + u(i,j+2) );
Lapu(i,j) = lapp/(c2*dt*dt);
ggux = 0.0; gguy = 0.0; ggvx = 0.0; ggvy = 0.0;
ggux = ggux + w(1)*(u(i+2,j+1)-u(i,j+1));
ggux = ggux + w(5)*(u(i+2,j+2)-u(i,j));
ggux = ggux + w(2)*0.0;
ggux = ggux + w(6)*(-1.0)*(u(i,j+2)-u(i+2,j));
ggux = ggux + w(3)*(-1.0)*(u(i,j+1)-u(i+2,j+1));
ggux = ggux + w(7)*(-1.0)*(u(i,j)-u(i+2,j+2));
ggux = ggux + w(4)*0.0;
ggux = ggux + w(8)*(u(i+2,j)-u(i,j+2));
ggvx = ggvx + w(1)*(v(i+2,j+1)-v(i,j+1));
ggvx = ggvx + w(5)*(v(i+2,j+2)-v(i,j));
ggvx = ggvx + w(2)*0.0;
ggvx = ggvx + w(6)*(-1.0)*(v(i,j+2)-v(i+2,j));
ggvx = ggvx + w(3)*(-1.0)*(v(i,j+1)-v(i+2,j+1));
ggvx = ggvx + w(7)*(-1.0)*(v(i,j)-v(i+2,j+2));
ggvx = ggvx + w(4)*0.0;
ggvx = ggvx + w(8)*(v(i+2,j)-v(i,j+2));
gguy = gguy + w(1)*0.0;
gguy = gguy + w(5)*(u(i+2,j+2)-u(i,j));
gguy = gguy + w(2)*(u(i+1,j+2)-u(i+1,j));
gguy = gguy + w(6)*(u(i,j+2)-u(i+2,j));
gguy = gguy + w(3)*0.0;
gguy = gguy + w(7)*(-1.0)*(u(i,j)-u(i+2,j+2));
gguy = gguy + w(4)*(-1.0)*(u(i+1,j)-u(i+1,j+2));
gguy = gguy + w(8)*(-1.0)*(u(i+2,j)-u(i,j+2));
ggvy = ggvy + w(1)*0.0;
ggvy = ggvy + w(5)*(v(i+2,j+2)-v(i,j));
ggvy = ggvy + w(2)*(v(i+1,j+2)-v(i+1,j));
ggvy = ggvy + w(6)*(v(i,j+2)-v(i+2,j));
ggvy = ggvy + w(3)*0.0;
ggvy = ggvy + w(7)*(-1.0)*(v(i,j)-v(i+2,j+2));
ggvy = ggvy + w(4)*(-1.0)*(v(i+1,j)-v(i+1,j+2));
ggvy = ggvy + w(8)*(-1.0)*(v(i+2,j)-v(i,j+2));
Gv_x(i,j) = ggvx/(2*c2*dt);
Gv_y(i,j) = ggvy/(2*c2*dt);
Gu_x(i,j) = ggux/(2*c2*dt);
Gu_y(i,j) = gguy/(2*c2*dt);
end
end
for i = 1:nx
for j = 1:ny
mu(i,j) = 4.0*beta*u(i+1,j+1)*(u(i+1,j+1)-1.0)*(u(i+1,j+1)+1.0) - kk*Lapu(i,j);
+ 1*( kk*v(i+1,j+1)*Lapu(i,j) + kk*( Gv_x(i,j)*Gu_x(i,j) + Gv_y(i,j)*Gu_y(i,j) ) + W*v(i+1,j+1)*u(i+1,j+1) );
if(v(i+1,j+1) >= 1.0-etaa)
mv(i,j) = kb*T*( log(v(i+1,j+1)) + 1.0 + (v(i+1,j+1)-1.0)/etaa + log(etaa) ) - beta*( u(i+1,j+1)^2-1 )^2 + 0.5*W*u(i+1,j+1)^2;
elseif( v(i+1,j+1) > etaa && v(i+1,j+1) < 1.0-etaa)
mv(i,j) = kb*T*( log( v(i+1,j+1) ) - log(1.0-v(i+1,j+1)) ) - beta*( u(i+1,j+1)^2-1 )^2 + 0.5*W*u(i+1,j+1)^2;
else
mv(i,j) = kb*T*( -log(1.0-v(i+1,j+1)) - 1.0 + v(i+1,j+1)/etaa +log(etaa) ) - beta*( u(i+1,j+1)^2-1 )^2 + 0.5*W*u(i+1,j+1)^2;
end
end
end
end
%stream.m
function [f] = stream(f)
f(:,:,1) = circshift( squeeze(f(:,:,1)), [+1,+0] );
f(:,:,2) = circshift( squeeze(f(:,:,2)), [+0,+1] );
f(:,:,3) = circshift( squeeze(f(:,:,3)), [-1,+0] );
f(:,:,4) = circshift( squeeze(f(:,:,4)), [+0,-1] );
f(:,:,5) = circshift( squeeze(f(:,:,5)), [+1,+1] );
f(:,:,6) = circshift( squeeze(f(:,:,6)), [-1,+1] );
f(:,:,7) = circshift( squeeze(f(:,:,7)), [-1,-1] );
f(:,:,8) = circshift( squeeze(f(:,:,8)), [+1,-1] );
end