LBM matlab.zip

% Oscillating % 入口密度边界，周期波动 clear;clc;close all; global cx cy w m n cx=[0,1,0,-1,0,1,-1,-1,1]; cy=[0,0,1,0,-1,1,1,-1,-1]; w=[4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36]; lx=301;ly=41; m=lx;n=ly; u0 = 0.1; % horizontal lid velocity v0 = 0; % vertical lid velocity Re = 50; % Reynolds number nu = u0 *ly / Re; % kinematic viscosity omega = 1. / (3*nu+1./2.); % relaxation parameter maxT = 30000; % total number of iterations [x,y]=meshgrid(1:m,1:n);[sx,sy]=meshgrid(1:16:m,1:16:n); u=u0*ones(m,n);v=zeros(m,n);rho=ones(m,n); fin=eqdf(u,v,rho); g=0.001*ones(m,n); tic for cycle=1:maxT f_eq=eqdf(u,v,rho); force=reshape(kron(w,rho).*kron(cx,g),m,n,9); % Collesion fout=omega*f_eq+(1-omega)*fin+force; % Streaming for k=1:9 fin(:,:,k)=circshift(fout(:,:,k),[cx(k),cy(k)]); end % BC % top 反弹边界上（2:m-1,n） fin(:,n,[ 5 8 9])=fin(:,n,[3 6 7]); % bottom 反弹边界 下（2:m-1,1） fin(:,1,[ 3 6 7])=fin(:,1,[5 8 9]); % left 压力入口边界 左（1,2:n-1） rho0=1.1+1*sin(2*pi*cycle/1000); ue=1-((sum(fin(1,2:n-1,[1 3 5]),3)+2*sum(fin(1,2:n-1,[4 7 8]),3))/rho0); fin(1,2:n-1,2)=fin(1,2:n-1,4)+2/3*rho0*ue; fin(1,2:n-1,6)=fin(1,2:n-1,8)-0.5*(fin(1,2:n-1,3)-fin(1,2:n-1,5))+rho0*ue/6; fin(1,2:n-1,9)=fin(1,2:n-1,7)+0.5*(fin(1,2:n-1,3)-fin(1,2:n-1,5))+rho0*ue/6; % Right 开放边界 右（m,2:n-1） fin(m,2:n-1,[ 4 7 8])=fin(m-1,2:n-1,[4 7 8]); % Macro Var rho=sum(fin,3); u=reshape(cx*(reshape(fin,[m*n 9]))',m,n)./rho; v=reshape(cy*(reshape(fin,[m*n 9]))',m,n)./rho; u(1,2:n-1)=u0; v(1,2:n-1)=0; % imagesc(sqrt(u.^2+v.^2)); if mod(cycle,10)==0 % subplot(2,2,1),plot(u(2,1:51));title('0.02');axis([1 51 -1 1]); % subplot(2,2,2),plot(u(32,1:51));title('0.3');axis([1 51 -1 1]); % subplot(2,2,3),plot(u(62,1:51));title('0.6');axis([1 51 -1 1]); % subplot(2,2,4),plot(u(2:201,25));title('0.9');axis([1 201 -1 1]); % % % plot(cycle/200,u(90,25),'ro');axis([1 51 -0.05 0.3]) meshc(u);axis([0 ,41,0,301,-1, 1]); new_level = -1; % Get the handle to each patch object h = findobj('type','patch'); % Create a loop to change the height of each contour zd = get(h,'ZData'); for i = 1:length(zd) set(h(i),'ZData',new_level*ones(length(zd{i}),1)) end view([113 26]); title(cycle/1000) %hold on % clf; % vv=cumsum(-0.5*(v(2:end,:)+v(1:end-1,:))); % uu=cumsum(0.5*(u(:,2:end)+u(:,1:end-1)),2); % uu(1,:)=[];vv(:,1)=[]; % strf=uu+vv; % strf=strf/max(max(abs(strf))); % contour(strf',10); % % %%% % axis equal off; drawnow toc end end %hold off