MATLAB code: Lattice Boltzmann method (D2Q9) for 2D Cahn-Hilliard equation

%Main code: Main.m

clear all;

nx = 101; ny = 101;

f = zeros(nx,ny,9); feq = zeros(nx,ny,9);

u = zeros(nx+2,ny+2); mu = zeros(nx,ny);

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 = 100.0; yl = 100.0;

dx = xl/(nx-1);

dy = yl/(ny-1);

dt = 1.0;

x = (0:dx:xl);

y = (0:dy:yl);

D = 3.0;

sigma = 0.1;

% kk = (3./8.)sigmaD;

% beta = (3./4.)*sigma/D;

beta = sqrt(0.09/(8*4.5));

kk = 4.5*beta;

eta = 1.0;

M = 1.0;

count = 0;

tau = 1.0; %M/(dtetac2)+0.5;

omega = 1./tau;

%setting initial condition

for i = 2:nx+1

for j = 2:ny+1

u(i,j) = -0.3 + 0.01*(2*rand-1);

end

end

figure(1);

result(nx,ny,x,y,u);

%Main Loop

while count < 4000

%Collitions

[f] = collition(nx,ny,u,mu,cx,cy,eta,omega,f,w);

%Streaming

[f] = stream(f);

%Boundary condition

[f] = boundary(nx,ny,f);

%Calculate u,mu

[u,mu] = ruv(nx,ny,u,mu,f,w,kk,beta,c2,dt);

count = count + 1

if(mod(count,400) == 0)

    figure(count);

    result(nx,ny,x,y,u);

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

%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

%result.m

function result(nx,ny,x,y,u)

xl = linspace(0,100,nx);yl = xl;

[xxl,yyl]=meshgrid(xl,yl);

A = u’;

surf(xxl,yyl,A(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] = ruv(nx,ny,u,mu,f,w,kk,beta,c2,dt)

Lapu = zeros(nx,ny);

u(2:nx+1,2:ny+1) = sum(f,3);

%periodic condition

u(1,:) = u(nx+1,:);

u(nx+2,:) = u(2,:);

u(:,1) = u(:,ny+1);

u(:,ny+2) = u(:,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);

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

 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