MATLAB code: Lattice Boltzmann method (D2Q9) for fluid flow with heat transfer

%Main code

clear all;

nx = 81; ny = 81;

f = zeros(nx,ny,9); g = zeros(nx,ny,9); rhog = ones(nx,ny);

u = zeros(nx,ny); v = zeros(nx,ny);

rho = ones(nx,ny); x = zeros(nx); y = zeros(ny);

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;

dx = 1.0; dy = 1.0;

xl = 1.0; yl = 1.0;

dx = xl/(nx-1);

dy = yl/(ny-1);

x = (0:dx:xl);

y = (0:dy:yl);

alpha = 0.05;

omega = 1.0/(3.0*alpha+0.5);

pr = 0.71;

Ra = 1.0e5;

alphag = alpha/pr;

ym1 = real(ny-1);

gbeta = Raalphaalphag/(ym1ym1ym1);

omegag = 1.0/(3.0*alphag+0.5);

count = 0;

%Main Loop

while count < 70000

%Collition

[f] = collition(nx,ny,u,v,cx,cy,omega,f,rho,w,rhog,gbeta);

%Streaming

[f] = stream(f);

%Boundary condition

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

%Calculate rho, u, v

[rho,u,v] = ruv(nx,ny,f);

%Temperature

[g] = gcol(nx,ny,u,v,cx,cy,omegag,g,rhog,w);

[g] = stream(g);

[g] = gbound(nx,ny,w,g);



rhog = sum(g,3);

count = count + 1;

end

result(nx,ny,x,y,u,v,rho,rhog);

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

%left boundary, bounce back

f(1,:,1) = f(1,:,3);

f(1,:,5) = f(1,:,7);

f(1,:,8) = f(1,:,6);

%right boundary, bounce back

f(nx,:,3) = f(nx,:,1);

f(nx,:,7) = f(nx,:,5);

f(nx,:,6) = f(nx,:,8);

%bottom boundary, bounce back

f(:,1,2) = f(:,1,4);

f(:,1,5) = f(:,1,7);

f(:,1,6) = f(:,1,8);

%top boundary, bounce back

f(:,ny,4) = f(:,ny,2);

f(:,ny,8) = f(:,ny,6);

f(:,ny,7) = f(:,ny,5);

end

function [f] = collition(nx,ny,u,v,cx,cy,omega,f,rho,w,rhog,gbeta)

for j = 1:ny

for i = 1:nx

    t1 = u(i,j)*u(i,j) + v(i,j)*v(i,j);

    for k = 1:9

        force = 3.0*w(k)*gbeta*rhog(i,j)*rho(i,j)*cy(k);

        if(i==1 || i==nx) force = 0.0;

        end

        if(j==1 || j==ny) force = 0.0;

        end

        t2 = u(i,j)*cx(k) + v(i,j)*cy(k);

        feq = rho(i,j)*w(k)*(1.0+3.0*t2+4.5*t2*t2-1.5*t1);

        f(i,j,k) = (1.0-omega)*f(i,j,k)+omega*feq+force;

    end

end

end

end

function [g] = gbound(nx,ny,w,g)

twall = 1.0;

%left boundary, T = 1.0

g(1,:,1) = twall*(w(1)+w(3))-g(1,:,3);

g(1,:,5) = twall*(w(5)+w(7))-g(1,:,7);

g(1,:,8) = twall*(w(8)+w(6))-g(1,:,6);

%right boundary, T = 0.0

g(nx,:,3) = -g(nx,:,1);

g(nx,:,7) = -g(nx,:,5);

g(nx,:,6) = -g(nx,:,8);

%bottom boundary, Insulated

g(:,1,2) = g(:,2,4);

g(:,1,5) = g(:,2,7);

g(:,1,6) = g(:,2,8);

%Top boundary, Insulated

g(:,ny,4) = g(:,ny-1,2);

g(:,ny,8) = g(:,ny-1,6);

g(:,ny,7) = g(:,ny-1,5);

end

function [g] = gcol(nx,ny,u,v,cx,cy,omegag,g,rhog,w)

for j = 1:ny

for i = 1:nx

 

    for k = 1:9

        t2 = u(i,j)*cx(k) + v(i,j)*cy(k);

        geq(i,j,k) = rhog(i,j)*w(k)*(1.0+3.0*t2);

        g(i,j,k) = (1.0-omegag)*g(i,j,k)+omegag*geq(i,j,k);

    end

end

end

end

function result(nx,ny,x,y,u,v,rho,rhog)

for j = 1:ny

uvm(j) = u((nx-1)/2,j);

tmj(j) = rhog((nx-1)/2,j);

end

for i = 1:nx

vvm(i) = v(i,(ny-1)/2);

end

figure

plot(tmj,y,’LineWidth’,2)

figure

plot(uvm,y,’LineWidth’,2)

xlabel(‘U’)

ylabel(‘Y’)

figure

plot(x,vvm,’LineWidth’,1.5)

xlabel(‘X’)

ylabel(‘V’)

figure

quiver(x,y,u,v,10)

%Stream function calculation

str = zeros(nx,ny);

for i = 1:nx

for j = 2:ny

    str(i,j) = str(i,j-1)+0.25*(rho(i,j)+rho(i,j-1))*(u(i,j)+u(i,j-1));

end

end

figure

contour(x,y,str);

end

function [rho,u,v] = ruv(nx,ny,f)

rho = sum(f,3);

for i = 1:nx

rho(i,ny) = f(i,ny,9)+f(i,ny,1)+f(i,ny,3)+2.0*(f(i,ny,2)+f(i,ny,6)+f(i,ny,5));

end

u = ( sum(f(:,:,[1 5 8]),3) - sum(f(:,:,[3 6 7]),3) )./rho;

v = ( sum(f(:,:,[2 5 6]),3) - sum(f(:,:,[4 7 8]),3) )./rho;

end

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