expandflow.rar

% expand flow simulation 2008-05-07 % % D2Q9模型 % C6 C2 C5 ^ j v n % \ | / | % C3 C9 C1 | % / | \ | % C7 C4 C8 |___________>i u m % function expandflow %----------------------------------------- % 规定初始计算变量 %----------------------------------------- clear all clc; tic; Re=200; %雷诺数？？ U=0.05; %上边界速度 m=150; n=50; %格子数 tau=0.5*(6*U*(n-1)/Re+1); %松弛时间 omega=1/tau %----------------------------------------- % 初始化密度 rho=ones(m,n); % 初始化速度(u,v) u=zeros(m,n); v=zeros(m,n); u(1,fix(n/4)+1:n-fix(n/4))=U; f=zeros(m,n,9); % 由已知条件给f赋初值 feq(:,:,1)=rho./9.0.*(1.0+3.0.*u+4.5.*u.^2-1.5.*(u.^2+v.^2)); feq(:,:,2)=rho./9.0.*(1.0+3.0.*v+4.5.*v.^2-1.5.*(u.^2+v.^2)); feq(:,:,3)=rho./9.0.*(1.0-3.0.*u+4.5.*u.^2-1.5.*(u.^2+v.^2)); feq(:,:,4)=rho./9.0.*(1.0-3.0.*v+4.5.*v.^2-1.5.*(u.^2+v.^2)); feq(:,:,5)=rho./36.0.*(1.0+3.0.*(u+v)+4.5.*(u+v).^2-1.5.*(u.^2+v.^2)); feq(:,:,6)=rho./36.0.*(1.0+3.0.*(-u+v)+4.5.*(-u+v).^2-1.5.*(u.^2+v.^2)); feq(:,:,7)=rho./36.0.*(1.0+3.0.*(-u-v)+4.5.*(-u-v).^2-1.5.*(u.^2+v.^2)); feq(:,:,8)=rho./36.0.*(1.0+3.0.*(u-v)+4.5.*(u-v).^2-1.5.*(u.^2+v.^2)); feq(:,:,9)=4.0.*rho./9.0.*(1.0-1.5.*(u.^2+v.^2)); f=feq; % 设置边界 bound=zeros(m,n); bound(fix(m/4):m,[1 n])=1; bound(1:fix(m/4),[1:fix(n/4) n-fix(n/4)+1:n])=1; balk=find(bound); size=m*n; step=0:size:7*size; to_reflect=[balk+step(1) balk+step(2) balk+step(3) balk+step(4) ... balk+step(5) balk+step(6) balk+step(7) balk+step(8)]; reflect=[balk+step(3) balk+step(4) balk+step(1) balk+step(2) ... balk+step(7) balk+step(8) balk+step(5) balk+step(6)]; avu=1; prevavu=1; t=0; numactivenodes=sum(sum(1-bound)); % 迭代计算 while (t<30000&1e-10<abs((prevavu-avu)/avu))|t<1000 %------------------------------------- % 1）输运； % 2）求解此时的宏观变量； % 3）由上步所得宏观变量求解该状态下的平衡分布函数； % 4）碰撞处理，即Boltzmann方程； % 5）边界处理； %------------------------------------ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 1、输运 f(m,:,:)=f(m-1,:,:); %下边界处完全发展 % z=1,(x-1,y)-- rel='nofollow' onclick='return false;'>(x,y) f(2:m,:,1)=f(1:m-1,:,1); % z=2,(x,y-1)-->(x,y) f(:,:,2)=f(:,[n 1:n-1],2); % z=3,(x+1,y)-->(x,y) f(1:m-1,:,3)=f(2:m,:,3); % z=4,(x,y+1)-->(x,y) f(:,:,4)=f(:,[2:n 1],4); % z=5,(x-1,y-1)-->(x,y) f(2:m,:,5)=f(1:m-1,[n 1:n-1],5); % z=6,(x+1,y-1)-->(x,y) f(1:m-1,:,6)=f(2:m,[n 1:n-1],6); % z=7,(x+1,y+1)-->(x,y) f(1:m-1,:,7)=f(2:m,[2:n 1],7); % z=8,(x-1,y+1)-->(x,y) f(2:m,:,8)=f(1:m-1,[2:n 1],8); % z=9,(x,y)-->(x,y)，即静止不动 % f(:,:,9)=f(:,:,9); %经输运后得到 f %------------------------------------ %反弹边界处理 boundback=f(to_reflect); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 2、求解宏观变量 %密度 rho=sum(f,3); %速度 u=(sum(f(:,:,[1 5 8]),3)-sum(f(:,:,[3 6 7]),3))./rho; v=(sum(f(:,:,[2 5 6]),3)-sum(f(:,:,[4 7 8]),3))./rho; %边界处限制条件 rho(balk)=0.000; u(balk)=0.000; v(balk)=0.000; rho(1,fix(n/4)+1:n-fix(n/4))=1.000; u(1,fix(n/4)+1:n-fix(n/4))=U; v(1,fix(n/4)+1:n-fix(n/4))=0.000; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 3、求解平衡分布函数 feq(:,:,1)=rho./9.0.*(1.0+3.0.*u+4.5.*u.^2-1.5.*(u.^2+v.^2)); feq(:,:,2)=rho./9.0.*(1.0+3.0.*v+4.5.*v.^2-1.5.*(u.^2+v.^2)); feq(:,:,3)=rho./9.0.*(1.0-3.0.*u+4.5.*u.^2-1.5.*(u.^2+v.^2)); feq(:,:,4)=rho./9.0.*(1.0-3.0.*v+4.5.*v.^2-1.5.*(u.^2+v.^2)); feq(:,:,5)=rho./36.0.*(1.0+3.0.*(u+v)+4.5.*(u+v).^2-1.5.*(u.^2+v.^2)); feq(:,:,6)=rho./36.0.*(1.0+3.0.*(-u+v)+4.5.*(-u+v).^2-1.5.*(u.^2+v.^2)); feq(:,:,7)=rho./36.0.*(1.0+3.0.*(-u-v)+4.5.*(-u-v).^2-1.5.*(u.^2+v.^2)); feq(:,:,8)=rho./36.0.*(1.0+3.0.*(u-v)+4.5.*(u-v).^2-1.5.*(u.^2+v.^2)); feq(:,:,9)=4.0.*rho./9.0.*(1.0-1.5.*(u.^2+v.^2)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 4、碰撞处理 f=(1-omega).*f+omega.*feq; %Boltzmann方程 f; %经碰撞处理后后所得f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 5、边界条件 %反弹边界条件 f(reflect)=boundback; f(find(f<0))=0.0000; % f>=0 prevavu=avu; avu=sum(sum(u))/numactivenodes; t=t+1; if rem(t,200)==0 t end end %密度 rho=sum(f,3); %速度 u=(sum(f(:,:,[1 5 8]),3)-sum(f(:,:,[3 6 7]),3))./rho; v=(sum(f(:,:,[2 5 6]),3)-sum(f(:,:,[4 7 8]),3))./rho; %界处限制条件 rho(balk)=0.000; u(balk)=0.000; v(balk)=0.000; rho(1,fix(n/4)+1:n-fix(n/4))=1.000; u(1,fix(n/4)+1:n-fix(n/4))=U; v(1,fix(n/4)+1:n-fix(n/4))=0.000; title='expand flow simulation'; SpeedDataOutput(title,u,v); figure;colormap(gray(2));image(2-bound');hold on; quiver(1:m,1:n,u',v'); axis equal; axis xy; toc