采用格子bolzmann方法MATLAB模拟圆柱绕流

function cylinder_LBM
% cylinder.m: Flow around a cyliner using LBM                 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 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      = 40000;      % total number of iterations 
tPlot    = 5;                % cycles 
 
% D2Q9 LATTICE CONSTANTS 
t    = [4/9 1/91/91/91/9 1/361/361/361/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 = [2:（ly-1）]; 
 
[yx] = meshgrid（1:ly1:lx）; 
obst = （x-obst_x）.^2 + （y-obst_y）.^2 <= obst_r.^2; 
obst（:[1ly]） = 1; 
bbRegion = find（obst）; 
 
% INITIAL CONDITION: （rho=0 u=0） ==> fIn（i） = t（i） 
fIn = reshape（ t‘ * ones（1lx*ly） 9 lx ly）; 
 
% MAIN LOOP （TIME CYCLES） 
for cycle = 1:max