Cavity problem by Lattice-Boltzmann method

Hi every body, I have written a MATLAB code for Lid-Driven cavity problem by Lattice-Boltzmann method. Although it seems that my code is OK, I could not get the correct figure. I spent one month or more and I got crazy! By the way, for seeing the streamlines, is it correct to use "contour" command (the last line in my code)?! I would greatly appriciate any help. my code is as follows:

if true
%%%%%%%%%%%%%%%
% Script file: LidDrivenCavity.m
%                                                                        
%  Lattice structure:
%  c4  c3   c2 
%    \  |  /          
%  c5 -c9 - c1 
%    /  |  \          
%  c6  c7   c8 
%

tic; hold on; clc; clear; nx=100; ny=100; tstep=40000; alpha=0.01; omega=1.0; u_ini=0.1; v_ini=0; Re=u_ini*nx/alpha w1=4/9; w2=1/9; w3=1/36; f=ones(nx,ny,9); f_eq=f; density=2.7;

for ii=1:tstep

    % Propegate (This part of code [propegate] is always constant for all LBM
    % problems.)
    f(:,:,4)=f([2:nx 1],[ny 1:ny-1],4); f(:,:,3)=f(:,[ny 1:ny-1],3);
    f(:,:,2)=f([nx 1:nx-1],[ny 1:ny-1],2); f(:,:,5)=f([2:nx 1],:,5);
    f(:,:,1)=f([nx 1:nx-1],:,1); f(:,:,6)=f([2:nx 1],[2:ny 1],6);
    f(:,:,7)=f(:,[2:ny 1],7); f(:,:,8)=f([nx 1:nx-1],[2:ny 1],8);

    % Boundary Conditions
    %At i=1, Bounceback
    f(1,:,1)=f(1,:,5); f(1,:,2)=f(1,:,6); f(1,:,8)=f(1,:,4);
    %At i=nx, Bounceback
    f(nx,:,4)=f(nx,:,8); f(nx,:,5)=f(nx,:,1); f(nx,:,6)=f(nx,:,2);
    %At j=1, Bounceback
    f(:,1,2)=f(:,1,6); f(:,1,3)=f(:,1,7); f(:,1,4)=f(:,1,8);
    %At j=ny, Know Velocity
    densityN=f(:,ny,9)+f(:,ny,1)+f(:,ny,5)+2*(f(:,ny,3)+f(:,ny,4)+f(:,ny,2));
    f(:,ny,7)=f(:,ny,3);
    f(:,ny,6)=f(:,ny,2)+0.5*(f(:,ny,1)-f(:,ny,5))-u_ini.*densityN/2;
    f(:,ny,8)=f(:,ny,4)+0.5*(f(:,ny,5)-f(:,ny,1))+u_ini.*densityN/2;

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

    f_eq(:,:,1)=w2*density.*(1+3*u+9/2*u.^2-3/2*(u.^2+v.^2));
    f_eq(:,:,3)=w2*density.*(1+3*v+9/2*v.^2-3/2*(u.^2+v.^2));
    f_eq(:,:,5)=w2*density.*(1-3*u+9/2*u.^2-3/2*(u.^2+v.^2));
    f_eq(:,:,7)=w2*density.*(1-3*v+9/2*v.^2-3/2*(u.^2+v.^2));

    f_eq(:,:,2)=w3*density.*(1+3*(u+v)+9/2*(u+v).^2-3/2*(u.^2+v.^2));
    f_eq(:,:,4)=w3*density.*(1+3*(-u+v)+9/2*(-u+v).^2-3/2*(u.^2+v.^2));
    f_eq(:,:,6)=w3*density.*(1-3*(u+v)+9/2*(u+v).^2-3/2*(u.^2+v.^2));
    f_eq(:,:,8)=w3*density.*(1+3*(u-v)+9/2*(u-v).^2-3/2*(u.^2+v.^2));

    f_eq(:,:,9)=w1*density.*(1-3/2*(u.^2+v.^2));

    f=omega*f_eq+(1-omega)*f;

end

contour(u',100); toc; end

Answer 1

[Aravind Baskar]

Although the streamlines are fine, error is not reducing after each time step, it is fluctuating.

[sthavishtha]

the fluctuation of errors, called as numerical instability is the main problem of SRT model of Lattice Boltzmann Method (method used here). In such cases, its recommended to adopt Multiple Relaxation Time (MRT) or Two-Relaxation time (TRT) models of Lattice Boltzmann Method.

Answer 2

[sthavishtha]

the method which you have suggested for plotting streamlines actually plots the u-velocity contours.For observing the streamlines, you can use streamslice, though its mostly preferred for 3D.

[x,y]=meshgrid(1:nx,1:ny); streamslice((x-1)/nx,(y-1)/ny,u',v');

Though the quality of the plotted streamlines is not great, you can also export the velocity data to a specific data file for visualization in widely used softwares - Tecplot, Visit and Paraview. Alternatively, if you want to use Matlab only, you may calculate the streamfunction directly from flowfun(u,v) - from the link attached. In that case, you may directly plot the contour of streamfunction.

Hope that helps.

Calculation of vorticity and stream function