LBM模拟二维平板流matlab代码

%    ------------------
%    ->                   二维通道，左进右出，上下为壁面
%    ------------------
% D2Q9模型
% C6  C2  C5    
%   \  |  /     
% C3  C9  C1    %原C0在此表示为C9
%   /  |  \     
% C7  C4  C8   

clear all
clc
m=41;
n=1001;
dt=1;
dx=1;
dy=1;
x=0:dx:1000;
y=0:dy:40;
alpha=0.02;%粘度
omega = 1/（alpha*3+0.5）;
uo=0.2;%初速度
rhoo=5;%初始密度
mstep=40000;

%初始化
rho=ones（mn）.*rhoo;
u=zeros（mn）;
v=zeros（mn）;
u（:1）=uo;
w=[1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36 4/9];  %原C0在此表示为C9
for i=1:9
    f（::i）=w（i）.*rho;
end
cx = [1 0 -1 0 1 -1 -1 1 0];
cy = [0 1 0 -1 1 1 -1 -1 0];
tic
for kk=1:mstep
    %碰撞过程
    t1 = u.^2+v.^2;
    for i = 1:9
        t2 = u.*cx（i） + v.*cy（i）;
        feq（::i）=w（i）.*rho.*（1+3*t2+4.5*t2.^2-1.5*t1）;
        f（::i）=（1-omega）.*f（::i）+omega.*feq（::i）;
    end
    
    %迁移过程
    f（:2:n1）=f（:1:n-11）;
    f（2:m:2）=f（1:m-1:2）;
    f（:1:n-13）=f（:2:n3）;
    f（1:m-1:4）=f（2:m:4）;
    f（2:m2:n5）=f（1:m-11:n-15）;
    f（2:m1:n-16）=f（1:m-12:n6）;
    f（1:m-11:n-17