Poiseuille_flow.zip

% Poiseuille_flow 直接作用力形式 压力驱动形式不了解 %目前看起来是对的 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=81;ly=21; m=lx;n=ly; u0 = 0; % horizontal lid velocity v0 = 0; % vertical lid velocity Re = 100; % Reynolds number nu = 0.1; % kinematic viscosity %% omega = 1. / (3*nu+1./2.); % relaxation parameter maxT = 30000; % total number of iterations u=u0*ones(m,n);v=zeros(m,n);rho=ones(m,n); fin=eqdf(u,v,rho); acc=0.001; F=zeros(m,n,9); tic for cycle=1:maxT f_eq=eqdf(u,v,rho); % Collesion fout=omega*f_eq+(1-omega)*fin; for i=1:9; %fout(:,:,i)=fout(:,:,i)+ 3*w(i)*(cx(i)*acc*rho); % LGA方式 Fi=((cx(i)-u)*acc+(cy(i)-v)*0).*f_eq(:,:,i)*(1-omega/2)*3; % He_Shan_Doolen方式 fout(:,2:n-1,i)=fout(:,2:n-1,i)+Fi(:,2:n-1); end % 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]); % Macro Var rho=sum(fin,3); u=reshape(cx*(reshape(fin,[m*n 9]))',m,n)./rho+0.5*acc; v=reshape(cy*(reshape(fin,[m*n 9]))',m,n)./rho; clf; % imagesc(sqrt(u.^2+v.^2)); if mod(cycle,50)==0 quiver(u',v',2); %streamline(stream2(x,y,u',v',sx,sy,[1 100])); %axis equal off; drawnow toc end end