• 迟子建
    了解作者
  • matlab
    开发工具
  • 1KB
    文件大小
  • rar
    文件格式
  • 0
    收藏次数
  • 1 积分
    下载积分
  • 3
    下载次数
  • 2018-11-12 03:07
    上传日期
扩管流 d2q9模型 matlab 距、
expandflow.rar
  • expandflow.m
    5.1KB
内容介绍
% 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
评论
    相关推荐
    • LBM_D2Q9_poiseuille.zip
      matlab 泊肃叶流动 LBM方法,很好用的小程序
    • LBM_D2Q9_RK_MODEL.zip
      lattice boltzmann 方法模拟两相流,采用RK模型
    • LBM_D2Q9_RK_MRT EE_MODEL.zip
      适合初学RK型格子玻尔兹曼模型,很好用,值得下载学习
    • (LBM)D2Q9.rar
      用格子玻尔兹曼方法实现二维流动传热模拟,移动边界方腔流。
    • d2q9.rar
      何老师教材上的LBM模拟顶盖驱动流的c++代码
    • D2Q9diffusion.zip
      二维正方形平板中的热扩散程序,使用lbm预测板中的等温线
    • LBM_D2Q9.zip
      LBM二维程序,基于D2Q9模型,用于计算方腔驱动流
    • D2Q9_2.rar
      模拟二维流体管柱绕流速度场,采用格子玻尔兹曼方法
    • d2q9.rar
      格子blotzmann方法的D2Q9模型,包括粒子的迁移和碰撞。
    • LBMVAL_D2Q9.rar
      应用LBM,采用D2Q9模型,用来模拟流动情况