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