testcouette.zip_LBM matlab_testcouette_格子玻尔兹曼

clear all clc m=41; n=1001; 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=5000; %初始化 rho=ones(m,n)*rhoo; u=zeros(m,n); v=zeros(m,n); 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 u_old=u; 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:n,1)=f(:,1:n-1,1);%right f(2:m,:,2)=f(1:m-1,:,2);%top f(:,1:n-1,3)=f(:,2:n,3);%left f(1:m-1,:,4)=f(2:m,:,4);%bottom f(2:m,2:n,5)=f(1:m-1,1:n-1,5); f(2:m,1:n-1,6)=f(1:m-1,2:n,6); f(1:m-1,1:n-1,7)=f(2:m,2:n,7); f(1:m-1,2:n,8)=f(2:m,1:n-1,8); %边界状态更新 rhon=f(m,:,9)+f(m,:,1)+f(m,:,3)+2*(f(m,:,2)+f(m,:,5)+f(m,:,6)); f(1,:,2)=f(1,:,4); f(1,:,5)=f(1,:,7); f(1,:,6)=f(1,:,8); f(m,:,4)=f(m,:,2); f(m,:,7)=f(m,:,5)-rhon*uo/6; f(m,:,8)=f(m,:,6)+rhon*uo/6; %f(:,1,1)=f(:,n,1); %f(:,1,5)=f(:,n,5); %f(:,1,8)=f(:,n,8); %f(:,n,3)=f(:,1,3); %f(:,n,6)=f(:,1,6); %f(:,n,7)=f(:,1,7); f(:,1,1)=f(:,2,1); f(:,1,5)=f(:,2,5); f(:,1,8)=f(:,2,8); f(:,n,3)=f(:,n-1,3); f(:,n,6)=f(:,n-1,6); f(:,n,7)=f(:,n-1,7); %计算宏观变量 rho = sum(f,3); rho(m,:)=f(m,:,9)+f(m,:,1)+f(m,:,3)+2*(f(m,:,2)+f(m,:,5)+f(m,:,6)); usum = zeros(m,n); vsum = zeros(m,n); for i=1:9 usum = usum + f(:,:,i).*cx(i); vsum = vsum + f(:,:,i).*cy(i); end u = usum./rho; v = vsum./rho; du=u-u_old; error=du./u; if error<=5*1e-9 break end end toc figure u1=u(:,100)/uo; plot(y,u1,'*-') hold on u2=u(:,200)/uo; plot(y,u2,'o-') u3=u(:,700)/uo; plot(y,u3,'--') hold off figure uv=sqrt(u.^2+v.^2); contourf(1:n,1:m,uv) colorbar figure quiver(u,v)