MRT_MRT-LBM_mrt_MRTLB_LBMMRT_LBM

%?????? MRT-LBM ??????????????? %??????????????????. %?? MRT-LBM ?????????? --???2017 %??????????QQ?? ?????????293267908??????????????????? %????? clc clear close all %?????? xLen=60; %??????? yLen=30; %??????? Re=10;%??? uMax=0.1;%?????? uAve=2/3*uMax; %???????????? nu=uAve*(yLen-1)/Re;%???? Cs=sqrt(1/3); %???? tau=1/2+3*nu; %???? omega=1/tau;%???? step=10000; %?????? y1=1; y2=yLen;%yLen???????? GG=uMax/((y1-y2)^2/4); f_u=GG*(2*nu);%????,??????????? f_v=0; %????,??????????? rhoo=1; %????? checkStep=100;%?????? saveStep=20; %?????? % filePath=uigetdir('*.*','D:\MyStudy\MATLAB\YuBrian');%????????????? VSSum=[];%??????????-? save Step ????? VSSum2=[]; %??????????-? check Step ???????????? %D2Q9 ???? w=[1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36 4/9 ]; %??????? cx=[ 1 0 -1 0 1 -1 -1 1 0];%??? x ???? cy=[ 0 1 0 -1 1 1 -1 -1 0];%??? y ???? %?? Gram-Schmidt???? M=[1 1 1 1 1 1 1 1 1;-4 -1 -1 -1 -1 2 2 2 2;4 -2 -2 -2 -2 1 1 1 1;... 0 1 0 -1 0 1 -1 -1 1;0 -2 0 2 0 1 -1 -1 1;0 0 1 0 -1 1 1 -1 -1;... 0 0 -2 0 2 1 1 -1 -1;0 1 -1 1 -1 0 0 0 0;0 0 0 0 0 1 -1 1 -1;]; %?????? 0 ??????matlab??9???????M???[-1 -1]?? M=circshift(M,[-1 -1]); %??????-???? % s1(9)=0; % s1(3)=0; % s1(5)=0; % s1(7)=omega; % s1(8)=s1(7); % % s1(4)=1.2; % s1(6)=s1(4); % % s1(1)=omega; % s1(2)=s1(1)-0.1; % 1 2 3 4 5 6 7 8 0 s1=[omega omega-0.1 0 1.2 0 1.2 omega omega 0]; s=diag(s1); %????-???? Minv=inv(M); % Minv_s=Minv*s; Minv_s=M\s; % sum(sum(Minv_s-s11)); %???????? y=1:yLen; u_analy=GG*(y-y1).*(-y+y2); rho=zeros(xLen,yLen); u=zeros(xLen,yLen); v=zeros(xLen,yLen); f=zeros(xLen,yLen,9); Fa=zeros(xLen,yLen,9); Fi=zeros(xLen,yLen,9); m_eq=zeros(9,1); % ????????????????=1  ?????????? 0???????????=?????? t1=u.^2+v.^2; rho(:,:)=rhoo; for i=1:xLen for j=1:yLen % rho(i,j)=rhoo; % u(i,j)=0; % v(i,j)=0; % t1=u(i,j)^2+v(i,j)^2;%?? U*V ??,t1 ?????????? for k=1:9 t2=u(i,j)*cx(k)+v(i,j)*cy(k); %U*C ?? f(i,j,k)=rho(i,j)*w(k)*(1+3*t2+4.5*t2^2-1.5*t1(i,j));%  ???????  ????=??????                         % Fa(k,i,j)=0; %???????????????,????? 0 end end end %??? tstart = tic;%???? [meshX,meshY]=meshgrid(1:xLen,yLen:-1:1);%??????? % figure % pc=pcolor(meshX,meshY,ones(yLen,xLen));%????? % shading interp %???? % axis equal % ylim([1 yLen]) for kk=1:step step2=kk; %??????? mm=max(max(rho)); %?????? %???????? if mm>100 break; end %1 ???? for i=1:xLen for j=1:yLen %   %????????? fmiddle=f(i,j,:); fvector=fmiddle(:); m=M*fvector; %??????-reference: mproving the Stability of the... %Multiple-Relaxation-Time Lattice Boltzmann Method by a... %viscosity Counteracting Approach m_eq(9,1)=rho(i,j); m_eq(1,1)=rho(i,j)*(-2+3*(u(i,j)^2+v(i,j)^2)); m_eq(2,1)=rho(i,j)*(1-3*(u(i,j)^2+v(i,j)^2)); m_eq(3,1)=rho(i,j)*u(i,j); m_eq(4,1)=-rho(i,j)*u(i,j); m_eq(5,1)=rho(i,j)*v(i,j); m_eq(6,1)=-rho(i,j)*v(i,j); m_eq(7,1)=rho(i,j)*(u(i,j)^2-v(i,j)^2); m_eq(8,1)=rho(i,j)*u(i,j)*v(i,j); %{ Reference: Timm-Kruger m_eq(9,1)=rho(i,j); m_eq(1,1)=rho(i,j)*(1-3*(u(i,j)^2+v(i,j)^2)); m_eq(2,1)=rho(i,j)*(9*u(i,j)^2*v(i,j)^2-3*(u(i,j)^2+v(i,j)^2)+1); m_eq(3,1)=rho(i,j)*u(i,j); m_eq(4,1)=rho(i,j)*(3*u(i,j)^3-u(i,j)); m_eq(5,1)=rho(i,j)*v(i,j); m_eq(6,1)=rho(i,j)*(3*v(i,j)^3-v(i,j)); m_eq(7,1)=rho(i,j)*(u(i,j)^2-v(i,j)^2); m_eq(8,1)=rho(i,j)*u(i,j)*v(i,j); %} %????-reference Timm Kruger: Eq. 10.42 Fi(i,j,9)=0; Fi(i,j,1)=6*(1-0.5*s1(1))*(u(i,j)*f_u+v(i,j)*f_v); Fi(i,j,2)=-6*(1-0.5*s1(2))*(u(i,j)*f_u+v(i,j)*f_v); Fi(i,j,3)=(1-0.5*s1(3))*f_u; Fi(i,j,4)=-(1-0.5*s1(4))*f_u; Fi(i,j,5)=(1-0.5*s1(5))*f_v; Fi(i,j,6)=-(1-0.5*s1(6))*f_v; Fi(i,j,7)=2*(1-0.5*s1(7))*(u(i,j)*f_u-v(i,j)*f_v); Fi(i,j,8)=(1-0.5*s1(8))*(u(i,j)*f_v+v(i,j)*f_u); %???????—????????????????? Fimiddle=Fi(i,j,:); Fivector=Fimiddle(:); m_temp=M\s*(m-m_eq)-M\Fivector; %M\Fivector=inv(M)*Fivector ????Timm book Eq. 10.7 or 10.13 for k=1:9 f(i,j,k)=f(i,j,k)-m_temp(k); end end end % 2 ???? f2=f;%?????????????f2 for i=1:9 f(:,:,i) = circshift(f2(:,:,i), [cx(i),cy(i),0]); end %??????-??????? u(:,1)=0; v(:,1)=0; u(:,yLen)=0; v(:,yLen)=0; for i=1:xLen %??? rho(i,1)=(f(i,1,9)+f(i,1,1)+f(i,1,3)+2*(f(i,1,4)+f(i,1,7)+f(i,1,8)))/(1-v(i,1)); f(i,1,2)=f(i,1,4)+2/3*rho(i,1)*v(i,1); f(i,1,5)=f(i,1,7)-1/2*(f(i,1,1)-f(i,1,3))+1/6*rho(i,1)*v(i,1)+1/2*rho(i,1)*u(i,1); f(i,1,6)=f(i,1,8)+1/2*(f(i,1,1)-f(i,1,3))+1/6*rho(i,1)*v(i,1)-1/2*rho(i,1)*u(i,1); %??? rho(i,yLen)=(f(i,yLen,9)+f(i,yLen,1)+f(i,yLen,3)+2*(f(i,yLen,2)+f(i,yLen,6)+f(i,yLen,5)))/(1+v(i,yLen)); f(i,yLen,4)=f(i,yLen,2)-2/3*rho(i,yLen)*v(i,yLen); f(i,yLen,7)=f(i,yLen,5)+1/2*(f(i,yLen,1)-f(i,yLen,3))-1/6*rho(i,yLen)*v(i,yLen)-1/2*rho(i,yLen)*u(i,yLen); f(i,yLen,8)=f(i,yLen,6)-1/2*(f(i,yLen,1)-f(i,yLen,3))-1/6*rho(i,yLen)*v(i,yLen)+1/2*rho(i,yLen)*u(i,yLen); end %????? %{ for i=1:xLen for j=1:yLen rho(i,j)=sum(f(i,j,3)); %???????? end end %} rho(:,:)=sum(f,3); for i=1:xLen for j=1:yLen uSum=0; vSum=0; for k=1:9 uSum=uSum+f(i,j,k)*cx(k); vSum=vSum+f(i,j,k)*cy(k); end u(i,j)=uSum/rho(i,j); %????? v(i,j)=vSum/rho(i,j); end end %? save Step ????????????      %{ if mod(kk,saveStep)==0 temp=0; for i=1:xLen for j=1:yLen temp=temp+sqrt(u(i,j)^2+v(i,j)^2);%?????? end end temp=sum(sum(sqrt(u.^2+v.^2))); VSSum2=[VSSum2;temp]; u2=rot90(u); v2=rot90(v); speed=sqrt(u2.^2+v2.^2); for i=1:xLen for j=1:yLen speed(j,i)=sqrt(u2(j,i)^2+v2(j,i)^2); %????? end end set(pc,'cdata',speed); drawnow saveas(gcf,['D:\MyStudy\MATLAB\YuBrian','\',num2str(kk/saveStep),'.tif']) picLength=kk/saveStep; %??????? %? check Step ?????? if mod(kk,checkStep)==0 VSSum=[VSSum;temp]; %???????????????*???????????? if length(VSSum)>=2 && abs(VSSum(end)-VSSum(end-1))/VSSum(end-1)<=(1e-6) fprintf('????:%d \n',kk); fprintf('????:%d\n',abs(VSSum(end)-VSSum(end-1))/(xLen*yLen)); break; end end end %} end runtime = toc(tstart); %?????? %????? u2=rot90(u); v2=rot90(v); %{ for i=1:xLen for j=1:yLen speed(j,i)=sqrt(u2(j,i)^2+v2(j,i)^2); end end %} speed=sqrt(u2.^2+v2.^2); %????? figure imagesc(speed); shading interp ylim([1 yLen]) figure surf(meshX,meshY,speed); shading interp %?????????? error=zeros(xLen,1); for i=1:xLen error(i)=(sqrt(sum((u(i,:)-u_analy).^2)))./sqrt(sum(u_analy.^2)); end L2=1/xLen*sum(error); fprintf('L2 ??:%d\n',L2); ydim=1:yLen; udim=u(1,:); udim2=u_analy'; figure plot(ydim,udim,'or') hold on plot(ydim,udim2) hold off xlabel('????','fontsize',16) ylabel('????','fontsize',16