{"id":3159,"date":"2023-12-16T20:42:15","date_gmt":"2023-12-16T12:42:15","guid":{"rendered":"http:\/\/81.70.49.155\/?p=3159"},"modified":"2023-12-16T20:45:32","modified_gmt":"2023-12-16T12:45:32","slug":"3%e5%9c%a8%e9%80%9a%e9%81%93%e6%b5%81%e5%9f%ba%e7%a1%80%e4%b8%8a%e6%94%b9%e6%88%90%e5%9c%86%e6%9f%b1%e7%bb%95%e6%b5%81matlab%e4%bb%a3%e7%a0%81","status":"publish","type":"post","link":"http:\/\/81.70.49.155\/?p=3159","title":{"rendered":"[LBM-3]\u5728\u901a\u9053\u6d41\u57fa\u7840\u4e0a\u6539\u6210\u5706\u67f1\u7ed5\u6d41MATLAB\u4ee3\u7801"},"content":{"rendered":"

\n\t
\n%\u7ea2\u8272<\/span>\u4ee3\u7801\u4e3a\u589e\u52a0\u4ee3\u7801\uff1b\u7eff\u5e95<\/span><\/strong>\u4e3a\u5c4f\u853d\u7684\u4ee3\u7801\n<\/p>\n

\n\t%\u521d\u59cb\u5316\n<\/p>\n

clc; clear; close; <\/span>
\n% lx=25; ly=9; <\/span>
\n% Re=1.25;<\/span>
\n% yu=0.1;<\/span>
\n% fx0=Re*8*yu^2\/(ly-1)^3; <\/span><\/p>\n

lx=400; ly=81; <\/p>\n

%\u6784\u5efa\u58c1\u9762\u53ca\u5706\u67f1\u533a\u57df<\/span>
\nBound=zeros(lx,ly);<\/span>
\nBound(:,[1,ly])=1;<\/span>
\nxc=lx\/4; yc=ly\/2; Rc=6;<\/span>
\nfor i=floor(xc)-7:floor(xc)+7<\/span>
\n for j=floor(yc)-7:floor(yc)+7<\/span>
\n if ((i-xc)^2+(j-yc)^2)^0.5<=Rc<\/span>
\n Bound(i,j)=1;<\/span>
\n end<\/span>
\n end<\/span>
\nend<\/span><\/p>\n

Bin=[]; BFin0=[]; BFin=[];<\/span>
\nopposite=[1,4,5,2,3,8,9,6,7];<\/span>
\nfor i=1:lx<\/span>
\n for j=1:ly<\/span>
\n if Bound(i,j)==1<\/span>
\n Bin=[Bin,i+lx*(j-1)];<\/span>
\n for k=1:9<\/span>
\n BFin=[BFin,(i+lx*(j-1))+(opposite(k)-1)*lx*ly];<\/span>
\n BFin0=[BFin0,(i+lx*(j-1))+(k-1)*lx*ly];<\/span>
\n end <\/span>
\n end<\/span>
\n end<\/span>
\nend<\/span><\/p>\n

Re=100;<\/span>
\nUin=0.1;<\/span>
\nyu=2*Rc*Uin\/Re;<\/span>
\nfx0=0; <\/span>
\nW=[4\/9,1\/9,1\/9,1\/9,1\/9,1\/36,1\/36,1\/36,1\/36];
\nex=[0, 1, 0, -1, 0 1,-1, -1, 1];
\ney=[0, 0, 1, 0, -1, 1, 1, -1,-1];
\ncs=1\/sqrt(3);<\/p>\n

rho=ones(lx,ly);
\nu=zeros(lx,ly);
\nv=zeros(lx,ly);
\nfx=zeros(lx,ly)+fx0;
\nfy=zeros(lx,ly);
\nM1=W(1)*rho; M2=W(2)*rho; M3=W(6)*rho;
\nF= cat(3,M1,M2,M2,M2,M2,M3,M3,M3,M3);
\nFeq=F;<\/p>\n

%\u8fc1\u79fb\u4e0b\u6807
\nD2x=[lx,1:lx-1]; D2y=[1:ly];
\nD3x=[1:lx]; D3y=[ly,1:ly-1];
\nD4x=[2:lx,1]; D4y=[1:ly];
\nD5x=[1:lx]; D5y=[2:ly,1];
\nD6x=[lx,1:lx-1]; D6y=[ly,1:ly-1];
\nD7x=[2:lx,1]; D7y=[ly,1:ly-1];
\nD8x=[2:lx,1]; D8y=[2:ly,1];
\nD9x=[lx,1:lx-1]; D9y=[2:ly,1];<\/p>\n

%% \u4e3b\u5faa\u73af===========================================
\nfor tStep=1:15000
\n% \u8fc1\u79fb
\nF0(:,:,1)=F(:,:,1);
\nF0(:,:,2)=F(D2x,D2y,2);
\nF0(:,:,3)=F(D3x,D3y,3);
\nF0(:,:,4)=F(D4x,D4y,4);
\nF0(:,:,5)=F(D5x,D5y,5);
\nF0(:,:,6)=F(D6x,D6y,6);
\nF0(:,:,7)=F(D7x,D7y,7);
\nF0(:,:,8)=F(D8x,D8y,8);
\nF0(:,:,9)=F(D9x,D9y,9);
\nF=F0;<\/p>\n

%\u8ba1\u7b97\u5b8f\u89c2\u91cf
\nrho=sum(F,3);
\nu=(sum(F(:,:,[2,6,9]),3)-sum(F(:,:,[4,7,8]),3)).\/rho+fx\/2.\/rho;
\nv=(sum(F(:,:,[3,6,7]),3)-sum(F(:,:,[5,8,9]),3)).\/rho+fy\/2.\/rho;<\/p>\n

%\u8bbe\u7f6e\u5b8f\u89c2\u8fb9\u754c
\nrho(:,1)=rho(:,2); %\u4e0b\u58c1\u9762
\nrho(:,ly)=rho(:,ly-1); %\u4e0a\u58c1\u9762
\nu(Bin)=0; v(Bin)=0; %\u6240\u6709\u56fa\u4f53<\/span><\/p>\n

rho(1,:)=rho(2,:); %\u5165\u53e3<\/span>
\nu(1,:)=Uin; v(1,:)=0; %<\/span><\/p>\n

rho(lx,:)=rho(lx-1,:); %\u51fa\u53e3<\/span>
\nu(lx,:)=u(lx-1,:); v(lx,:)=v(lx-1,:);<\/span><\/p>\n

%\u8ba1\u7b97\u5e73\u8861\u6001\u51fd\u6570
\nfor i=1:9
\n M1=(ex(i)*u+ey(i)*v)\/cs^2;
\n M2=M1.^2\/2;
\n M3=(u.^2+v.^2)\/2\/cs^2;
\n Feq(:,:,i)=rho.*W(i).*(1+M1+M2-M3);
\nend<\/p>\n

%\u78b0\u649e
\nfor i=1:9
\n M1=3*(ex(i)-u).*fx+(ey(i)-v).*fy;
\n M2=9*(ex(i)*u+ey(i)*v).*(ex(i).*fx+ey(i).*fy);
\n force(:,:,i)=W(i)*(1-0.5\/(0.5+3*yu))*(M1+M2);
\nend
\n F=F-1\/(0.5+3*yu)*(F-Feq)+force;
\n
\n%\u5fae\u89c2\u8fb9\u754c\u8bbe\u7f6e
\nF(BFin0)=F(BFin); %\u534a\u6b65\u53cd\u5f39\u64cd\u4f5c<\/span>
\nF(1,:,:)=Feq(1,:,:)+(F(2,:,:)-Feq(2,:,:));<\/span>
\nF(lx,:,:)=Feq(lx,:,:)+(F(lx-1,:,:)-Feq(lx-1,:,:));<\/span>
\n% F(:,[1,ly],:)=F(:,[1,ly],opposite);<\/span>
\n% F(:,1,:)=Feq(:,1,:)+(F(:,2,:)-Feq(:,2,:));
\n% F(:,ly,:)=Feq(:,ly,:)+(F(:,ly-1,:)-Feq(:,ly-1,:));<\/p>\n

%\u540e\u5904\u7406
\nif mod(tStep,100)==0
\nclc;clf
\ntStep
\nU=(u.^2+v.^2).^0.5;
\nsubplot(1,2,1),imagesc(U');
\ncolormap jet
\naxis equal
\ncy=(ly-1)\/2;
\ny0=fx0\/2\/yu*(cy.^2-([0:ly-1]-cy).^2);
\nNi=floor(lx\/2);
\nsubplot(1,2,2),hold on,plot(y0,'k'),plot(u(Ni,:),'r'); hold off
\ndrawnow
\nend<\/p>\n

end % \u4e3b\u5faa\u73af\u7ed3\u675f============================<\/p>\n","protected":false},"excerpt":{"rendered":"