Cavity problem by Lattice-Boltzmann method

12 次查看（过去 30 天）
Emad
Emad 2014-1-29
评论： sthavishtha 2017-1-26
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

请先登录，再回答此问题。

回答（2 个）

Aravind Baskar
Aravind Baskar 2016-4-16
Although the streamlines are fine, error is not reducing after each time step, it is fluctuating.
  1 个评论
sthavishtha
sthavishtha 2017-1-26
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.

请先登录，再进行评论。

sthavishtha
sthavishtha 2017-1-26
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

类别

在 Help Center 和 File Exchange 中查找有关 Fluid Dynamics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by