【资源】基于LBM，绕流问题的Matlab源代码

ywang

我从原来论坛copy过来的，呵呵~~！


the Program for Channel flow past a cylinderical obstacle, using a LB method

下面是MATLAB程序,程序为网络所得,版权归程序所有者 Jonas Latt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cylinder.m: Channel flow past a cylinderical   
%         obstacle, using a LB method        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lattice Boltzmann sample, written in C
% Copyright (C) 2006 Jonas Latt
% Address: Rue General Dufour 24, 1211 Geneva 4, Switzerland
% E-mail: Jonas.Latt@cui.unige.ch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
% You should have received a copy of the GNU General Public
% License along with this program; if not, write to the Free
% Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
% Boston, MA 02110-1301, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

% GENERAL FLOW CONSTANTS
lx       = 250;
ly       = 51;
obst_x = lx/5+1;   % position of the cylinder; (exact
obst_y = ly/2+1;   % y-symmetry is avoided)
obst_r = ly/10+1; % radius of the cylinder
uMax = 0.02;     % maximum velocity of Poiseuille inflow
Re   = 100;     % Reynolds number
nu   = uMax * 2.*obst_r / Re;   % kinematic viscosity
omega = 1. / (3*nu+1./2.);     % relaxation parameter
maxT   = 400000;   % total number of iterations
tPlot = 5;     % cycles

% D2Q9 LATTICE CONSTANTS
t = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cx = [ 0,   1, 0, -1, 0,   1, -1, -1,   1];
cy = [ 0,   0, 1, 0, -1,   1,   1, -1, -1];
opp = [ 1,   4, 5, 2, 3,   8,   9,   6,   7];
col = [2ly-1)];

[y,x] = meshgrid(1:ly,1:lx);
obst = (x-obst_x).^2 + (y-obst_y).^2 <= obst_r.^2;
obst(:,[1,ly]) = 1;
bbRegion = find(obst);

% INITIAL CONDITION: (rho=0, u=0) ==> fIn(i) = t(i)
fIn = reshape( t' * ones(1,lx*ly), 9, lx, ly);

% MAIN LOOP (TIME CYCLES)
for cycle = 1:maxT

  % MACROSCOPIC VARIABLES
  rho = sum(fIn);
  ux = reshape ( ...
      (cx * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;
  uy = reshape ( ...
      (cy * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;
   
  % MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS
    % Inlet: Poiseuille profile
  L = ly-2; y = col-1.5;
  ux(:,1,col) = 4 * uMax / (L*L) * (y.*L-y.*y);
  uy(:,1,col) = 0;
  rho(:,1,col) = 1 ./ (1-ux(:,1,col)) .* ( ...
    sum(fIn([1,3,5],1,col)) + ...
    2*sum(fIn([4,7,8],1,col)) );
    % Outlet: Zero gradient on rho/ux
  rho(:,lx,col) = 4/3*rho(:,lx-1,col) - ...
            1/3*rho(:,lx-2,col);
  uy(:,lx,col) = 0;
  ux(:,lx,col) = 4/3*ux(:,lx-1,col) - ...
            1/3*ux(:,lx-2,col);

  % COLLISION STEP
  for i=1:9
    cu = 3*(cx(i)*ux+cy(i)*uy);
    fEq(i,:, = rho .* t(i) .* ...
      ( 1 + cu + 1/2*(cu.*cu) ...
          - 3/2*(ux.^2+uy.^2) );
    fOut(i,:, = fIn(i,:, - ...
      omega .* (fIn(i,:,:)-fEq(i,:,:));
  end

  % MICROSCOPIC BOUNDARY CONDITIONS
  for i=1:9
      % Left boundary
      fOut(i,1,col) = fEq(i,1,col) + ...
      18*t(i)*cx(i)*cy(i)* ( fIn(8,1,col) - ...
      fIn(7,1,col)-fEq(8,1,col)+fEq(7,1,col) );
      % Right boundary
      fOut(i,lx,col) = fEq(i,lx,col) + ...
      18*t(i)*cx(i)*cy(i)* ( fIn(6,lx,col) - ...
      fIn(9,lx,col)-fEq(6,lx,col)+fEq(9,lx,col) );
      % Bounce back region
      fOut(i,bbRegion) = fIn(opp(i),bbRegion);
  end

  % STREAMING STEP
  for i=1:9
    fIn(i,:,:) = ...
      circshift(fOut(i,:,:), [0,cx(i),cy(i)]);
  end

  % VISUALIZATION
  if (mod(cycle,tPlot)==0)
    u = reshape(sqrt(ux.^2+uy.^2),lx,ly);
    u(bbRegion) = nan;
    imagesc(u');
    axis equal off; drawnow
  end
end

[ 本帖最后由 ywang 于 2009-6-22 15:44 编辑 ]

haixing

沙发！里面有一些字变成“ ”了，呵呵！

yushulinfeng

为什么有那么多笑脸呢

ywang

发帖时候，插入表情，就能看到表情对应的是啥符号了，O(∩_∩)O

rjqu

heisenbergg

这个如果弄成附件就比较好了

ywang

原帖由 heisenbergg 于 2010-1-29 13:37 发表

                               
登录/注册后可看大图

这个如果弄成附件就比较好了

恩，对，有道理！

onesupeng

其实，大家可以调试调试，然后加点注解共享，对初学者更有用

onesupeng

老鸟好像不大感冒

zwn

原帖由 onesupeng 于 2010-3-1 07:34 发表

                               
登录/注册后可看大图

其实，大家可以调试调试，然后加点注解共享，对初学者更有用

支持

juguang

谢谢楼主，现在正需要

way8210

arlu

正巧再看，为什么进口速度初始化的时候 y = col-1.5;？

kavine1986

kavine1986

:loveliness: :)