方腔顶盖驱动流流场数值预测
方腔顶盖驱动流动
一、二、问题描述方腔顶盖驱动流动如图1所示的一个简化两维方腔（高，宽都等于L)，内部充满水分。
上表面为移动墙，非维化速度为ｕ／ｕ0 =1。
其他三面为固定墙。
试求方腔内水分流动状态。
u=1, v=0u=0, v=0u=0,v=0u=0, v=0图1常微分方程理论只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法.二、离散格式数值解法:求解所有的常微分方程 计算解函数 y(x) 在一系列节点a = x 0< x 1<…<x n = b 处的近似值),...,1()(n i x y y i i =≈节点间距为步长，通常采用等距节点，即取 hi = h (常数)。
步进式：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。
因此只需建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。
欧拉方法1(,) 0,1,...n n n n y y h f x y n +=+=几何意义在假设 y n = y (x n )，即第 n 步计算是精确的前提下，考虑公式或方法本身带来的误差: R n = y (x n +1)y n +1 , 称为局部截断误差.显式欧拉公式一阶向前差商近似一阶导数223111232()[()()()()][ (,)] ()()h n n n n n n n n n h n R y x y y x hy x y x O h y hf x y y x O h +++'''=-=+++-+''=+推导如下：隐式欧拉公式x n +1点向后差商近似导数 推导如下：1()()()n n n y x y x y x h+-'≈111()()() ()()(,)n n n n nn n n n n y x y x hy x y x y y x y y h f x y +++'≈+↑≈≈=+11()()()n n n y x y x y x h++-'≈11()()()()n n n n ny x y x hy x y x y ++'≈+↑≈几何意义设已知曲线上一点 P n (x n , y n ),过该点作弦线，斜率为(x n +1 , y n +1 ) 点的方向场f (x ,y )方向,若步长h 充分小，可用弦线和垂线x =x n +1的交点近似曲线与垂线的交点。
方腔环流的流场计算
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
院系：工学院 姓名：刘广 学号：11309018 日期：2014 年 05 月 23 号
方腔环流示意图
我们取网格做如下说明，取正方形网格，网格数量由用户输入，但是建议不超过 200，因 为超过 200 计算量会急剧增大，雷诺数也是由用户输入，建议不超过 5000，因为超过 5000 流 场已经呈现出非线性，开始出现偏差，雷诺数超过 10000 甚至根本不能算出结果，这是因为强 非线性情况下描述流场行为的方程已经不能做如上简化， 同样的， 在计算过程中我们采用超松 弛迭代法，下面对离散形式做一下说明： 我们用二阶精度的差商代替微商做以下说明，得
姓名：刘广
学号：11309018
日期：2014 年 05 月 23 号
程序执行时，我输入的是10E-3作为计算精度，对于程序的说明程序注释已经说的很清楚， 这里不再赘述。最终会在目录下生成OUTPUT.txt 。 四：实验结果 最终会在目录下生成OUTPUT.txt 。然后我们可以用Matlab软件对计算出来的数据进一步整 理，画成云图和画出流线，我们做一下说明。
Re n n [( in, j 1 in, j 1 )( in in1, j in1, j )( in 1, j i 1, j ) ( , j 1 i , j 1 )]} 4 (1 ) in ,j
边界值的点需要特殊处理，然后由内点差分格式顺序扫描，一次次迭代直到满足用户输入的 精度为止。边界条件如下：
( i 1, j 2 i , j i 1, j ) / h 2 ( i , j 1 2 i , j i , j 1 ) / h 2 i , j ( i , j 1 i , j 1 )(
格子Boltzmann方法三种边界格式的对比分析
格子Boltzmann方法三种边界格式的对比分析刘连国;杨帆;王宏光【摘要】采用格子Boltzmann方法(Lattice Boltzmann Method- LBM)对二维顶盖驱动方腔流动进行数值模拟.在计算中分别使用半步长反弹、非平衡反弹、以及非平衡外推三种边界处理格式,并得到了不同格式对应的流线分布,流函数最小值、涡心坐标、几何中心线速度分布等.通过将所得结果与基准解进行比较,就三种边界格式的计算效率,计算精度、以及计算稳定性等方面进行了讨论和分析,为LBM计算中边界格式的选择提供了有益的参考.【期刊名称】《机械研究与应用》【年(卷),期】2012(000)001【总页数】5页(P18-22)【关键词】格子Boltzmann方法;边界处理格式;半步长反弹格式;非平衡反弹格式;非平衡外推格式【作者】刘连国;杨帆;王宏光【作者单位】上海理工大学能源与动力工程学院,上海200093;上海理工大学能源与动力工程学院,上海200093;上海理工大学能源与动力工程学院,上海200093【正文语种】中文【中图分类】O357.11 引言格子Boltzmann方法(LBM)是近年来迅速发展的一种新型数值计算方法。
边界条件的处理是LBM实施中一项非常关键的内容。
实际计算表明:选取不同的边界条件会对数值计算的精度、稳定性以及效率产生很大影响。
作为LBM的一个基本问题，边界条件的处理一直是流体力学一个重要的研究方面。
根据边界条件的类型，可将之分为两类:压力边界和速度边界［1］，其中的速度边界又可细分为:平直边界和曲面边界。
笔者从经典的流体力学问题二维顶盖方腔流模拟入手，对三种平直边界格式进行对比和分析，为LBM计算中边界格式的选择提供了有益的参考。
2 二维九点格子Boltzmann模型目前最常用的格子Boltzmann模型为LBGK模型，通过引入“单一弛豫时间”来简化Boltzmann方程中碰撞项的计算［2］。
九点格子LBGK模型的演化方程为:式中:(x，t)是在t时刻、x处的平衡态分布函数;τ为单一弛豫时间因子;eα为网格点各方向上的粒子速度。
顶盖驱动流计算程序
!此?程ì序ò用?QUICK格?式?求ó解a顶￥盖?驱y动ˉ流ⅰ?!对?压1力求ó解a利?用?人?工¤压1缩?法ぁ?求ó解aprogram QUICK_cavityimplicit none!mx和ímy分?别纄表括?示?网?格?节ú点?数簓integer :: i,j,mx,my,numcharacter(len=15) :: name,name1!c2表括?示?人?工¤压1缩?系μ数簓的?平?方?,dt非?稳è态?前°后ó时骸?间?间?隔?!dx和ídy表括?示?网?格?间?距à,x1和íy1表括?示?边?长¤,err表括?示?判D断?人?工¤压1缩?法ぁ?求ó解a收?敛?的?标括?准?!u0表括?示?顶￥盖?运?动ˉ的?速ù度è,relax表括?示?人?工¤压1缩?系μ数簓real(kind=8) :: c2,relax,dt,dx,dy,x1,y1,err,abc,u0!density表括?示?流ⅰ?体?密ü度è,Viscosity表括?示?流ⅰ?体?粘3度è,Re表括?示?雷ぁ?诺μ数簓real(kind=8) :: density,Viscosity,Re!表括?示?t时骸?刻ì的?速ù度è和í压1力值μreal(kind=8),allocatable :: u(:,:),v(:,:),p(:,:)!表括?示?t+dt时骸?刻ì的?速ù度è和í压1力值μreal(kind=8),allocatable :: un(:,:),vn(:,:),pn(:,:)!最?终?各÷节ú点?上?的?速ù度è和í压1力值μreal(kind=8),allocatable :: uc(:,:),vc(:,:),pc(:,:)!定¨义?涡D量?和í流ⅰ?函ˉ数簓real(kind=8),allocatable:: flow_function_u(:,:),flow_function_v(:,:),vortex_function(:,:)!给?各÷节ú点?定¨义?坐?标括?矩?阵óreal(kind=8),allocatable :: x(:),y(:)open(3,file="parameter.dat")read(3,*) mx,my,x1,y1 !读á入?网?格?节ú点?数簓和í边?长¤read(3,*) relax,dt,u0 !读á入?压1缩?系μ数簓,时骸?间?间?隔?和í顶￥盖?移?动ˉ速ù度èread(3,*) density,Viscosity !读á入?密ü度è和í粘3度èclose(3)allocate(u(mx,my+1),v(mx+1,my),p(mx+1,my+1))allocate(un(mx,my+1),vn(mx+1,my),pn(mx+1,my+1))allocate(uc(mx,my),vc(mx,my),pc(mx,my))allocate(flow_function_u(mx,my),flow_function_v(mx,my),vortex_function(mx,my))allocate(x(mx),y(my))dx=x1/(mx-1)dy=y1/(my-1)num=0err=100.0c2=relax**2Re=u0*density*x1/Viscosity!对?流ⅰ?场?进?行D初?始?化ˉdo i=1,mx+1,1do j=1,my+1,1p(i,j)=1.0end doend dodo i=1,mx,1do j=1,my+1,1u(i,j)=0.0if(j==my+1) u(i,j)=3.0*u0/2.0if(j==my) u(i,j)=u0/2.0end doend dodo i=1,mx+1,1do j=1,my,1v(i,j)=0.0end doend do!利?用?人?工¤压1缩?法ぁ?求ó解a流ⅰ?场?值μdo while(err>0.00001.and.num<1000000)err=0.0!调獭?用?quick子哩?程ì序ò中D的?QUICK离?散Α?格?式?计?算?un和ívn的?值μcall quick(u0,u,v,p,mx,my,dx,dy,dt,density,Viscosity,un,vn)!调獭?用?calp子哩?程ì序ò求ó解a压1强?pn值μcall calp(un,vn,pn,mx,my,dx,dy,dt,c2,p)!校￡验é流ⅰ?场?信?息￠，?得?到?收?敛?判D断?准?则òerr，?同?时骸?更ü新?u、￠v、￠p !利?用?单蹋?位?时骸?间?间?隔?前°后ó时骸?刻ì之?间?的?差?值μ的?绝?对?值μ作痢?为a 判D据Ydo i=1,mx,1do j=1,my+1,1abc=abs(un(i,j)-u(i,j))/dtif(abc>err) err=abcu(i,j)=un(i,j)end doend dodo i=1,mx+1,1do j=1,my,1abc=abs(vn(i,j)-v(i,j))/dtif(abc>err) err=abcv(i,j)=vn(i,j)end doend dodo j=1,my+1,1abc=(abs(pn(i,j)-p(i,j))/c2)/dtif(abc>err) err=abcp(i,j)=pn(i,j)end doend dowrite(*,*) "error=",errnum=num+1write(*,*) numend do!---------循-环·结á束?-------------------------------------------------------do i=1,mx,1x(i)=(i-1)*dxend dodo j=1,my,1y(j)=(j-1)*dyend do!计?算?节ú点?上?的?涡D量?do i=2,mx-1,1do j=2,my-1,1vortex_function(i,j)=(un(i,j+1)-un(i,j))/dy-(vn(i+1,j)-vn(i,j))/dx end doend dodo j=1,my,1vortex_function(1,j)=(un(1,j+1)-un(1,j))/dy-(vn(2,j)-vn(1,j))/dxvortex_function(mx,j)=(un(mx,j+1)-un(mx,j))/dy-(vn(mx+1,j)-vn(mx,j))/dx end dodo i=2,mx-1,1vortex_function(i,1)=(un(i,2)-un(i,1))/dy-(vn(i+1,1)-vn(i,1))/dxvortex_function(i,my)=(un(i,my+1)-un(i,my))/dy-(vn(i+1,my)-vn(i,my))/dx end do!计?算?节ú点?上?的?速ù度è值μdo i=1,mx,1do j=1,my,1uc(i,j)=0.5*(u(i,j)+u(i,j+1))vc(i,j)=0.5*(v(i,j)+v(i+1,j))pc(i,j)=0.25*(p(i,j)+p(i,j+1)+p(i+1,j)+p(i+1,j+1))end doend do!计?算?节ú点?上?的?流ⅰ?函ˉ数簓flow_function_u=0.0flow_function_v=0.0do i=2,mx-1,1flow_function_u(i,j)=dy*uc(i,j)+flow_function_u(i,j-1)flow_function_v(i,j)=-dx*vc(i,j)+flow_function_u(i-1,j)end doend do!输?出?数簓据Y到?文?件t夹D中Dwrite(name,"(f10.4)") Reopen(10,file='Reout'//trim(adjustl(name))//'.dat')write(10,*) 'TITLE = "result"'write(10,*) 'VARIABLES = "X","Y","U","V","P"'write(10,*) 'ZONE I=',mx,'J=',my,'F=POINT'write(10,"(5(f15.8,5x))") ((x(i),y(j),uc(i,j),vc(i,j),pc(i,j),i=1,mx),j=1,my)close(10)write(name1,"(f10.4)") Reopen(20,file='Refunction'//trim(adjustl(name1))//'.dat')write(20,*) 'TITLE = "result"'write(20,*) 'VARIABLES = "X","Y","FLOW_FUNCTION_U","FLOW_FUNCTION_V","VORTEX_FUNCTION"'write(20,*) 'ZONE I=',mx,'J=',my,'F=POINT'write(20,"(5(f15.8,5x))")((x(i),y(j),flow_function_u(i,j),flow_function_v(i,j),vortex_function(i,j),i=1,mx),j=1, my)close(20)stopend!利?用?一?阶×迎?风?格?式?和íQUICK格?式?求ó解a速ù度è值μsubroutine quick(u0,u,v,p,mx,my,dx,dy,dt,density,Viscosity,un,vn)implicit noneinteger :: i,j,mx,myreal(kind=8) :: dx,dy,dt,Viscosity,density,u0real(kind=8) :: u(mx,my+1),v(mx+1,my),p(mx+1,my+1)real(kind=8) :: un(mx,my+1),vn(mx+1,my)real(kind=8) :: fw,fe,fs,fn,df,dw,de,ds,dnreal(kind=8) :: aw,aww,ae,aee,as,ass,an,ann,apreal(kind=8),external :: alfa!求ó解ax方?向ò的?速ù度èundo i=3,mx-2,1do j=3,my-1,1!求ó解aQUICK格?式?中D的?系μ数簓fw=0.5*density*(u(i-1,j)+u(i,j))*dyfe=0.5*density*(u(i,j)+u(i+1,j))*dyfs=0.5*density*(v(i,j-1)+v(i+1,j-1))*dxfn=0.5*density*(v(i,j)+v(i+1,j))*dxdf=fe-fw+fn-fsdw=Viscosity*dy/dxde=Viscosity*dy/dxds=Viscosity*dx/dydn=Viscosity*dx/dyaw=dw+0.75*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw))*fwaww=-0.125*alfa(fw)*fwae=de-0.375*alfa(fe)*fe-0.75*(1.0-alfa(fe))*fe-0.125*(1.0-alfa(fw))*fwaee=0.125*(1.0-alfa(fe))*feas=ds+0.75*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs))*fsass=-0.125*alfa(fs)*fsan=dn-0.375*alfa(fn)*fn-0.75*(1.0-alfa(fn))*fn-0.125*(1.0-alfa(fs))*fsann=0.125*(1.0-alfa(fn))*fnap=aw+aww+ae+aee+as+ass+an+ann+dfun(i,j)=u(i,j)+dt/(dx*dy)*(-ap*u(i,j)+aw*u(i-1,j)+ae*u(i+1,j)+as*u(i,j-1)+an*u(i,j+1)+& aww*u(i-2,j)+aee*u(i+2,j)+ass*u(i,j-2)+ann*u(i,j+2))-dt*(p(i+1,j)-p(i,j))/dx end doend do!----------------------------------------------------------------------!内ú层?边?界?由?一?阶×迎?风?离?散Α?格?式?计?算?得?到?j=2do i=3,mx-2,1call upbound_u(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,un)end doj=mydo i=3,mx-2,1call upbound_u(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,un)end doi=2do j=2,my,1call upbound_u(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,un)end doi=mx-1do j=2,my,1call upbound_u(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,un)end do!--------------------------------------------------------------------!计?算?外猘层?边?界?do i=2,mx-1,1un(i,1)=-un(i,2)un(i,my+1)=2.0*u0-un(i,my)end dodo j=1,my+1,1un(1,j)=0.0un(mx,j)=0.0end do!--------------------------------------------------------------------!求ó解ay方?向ò的?速ù度èvndo i=3,mx-1,1do j=3,my-2,1!求ó解aQUICK格?式?中D的?系μ数簓fw=0.5*density*(u(i-1,j)+u(i-1,j+1))*dyfe=0.5*density*(u(i,j)+u(i,j+1))*dyfs=0.5*density*(v(i,j-1)+v(i,j))*dxfn=0.5*density*(v(i,j)+v(i,j+1))*dxdf=fe-fw+fn-fsdw=Viscosity*dy/dxde=Viscosity*dy/dxds=Viscosity*dx/dydn=Viscosity*dx/dyaw=dw+0.75*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw))*fwaww=-0.125*alfa(fw)*fwae=de-0.375*alfa(fe)*fe-0.75*(1.0-alfa(fe))*fe-0.125*(1.0-alfa(fw))*fwaee=0.125*(1.0-alfa(fe))*feas=ds+0.75*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs))*fsass=-0.125*alfa(fs)*fsan=dn-0.375*alfa(fn)*fn-0.75*(1.0-alfa(fn))*fn-0.125*(1.0-alfa(fs))*fsann=0.125*(1.0-alfa(fn))*fnap=aw+aww+ae+aee+as+ass+an+ann+dfvn(i,j)=v(i,j)+dt/(dx*dy)*(-ap*v(i,j)+aw*v(i-1,j)+ae*v(i+1,j)+as*v(i,j-1)+an*v(i,j+1)+& aww*v(i-2,j)+aee*v(i+2,j)+ass*v(i,j-2)+ann*v(i,j+2))-dt*(p(i,j+1)-p(i,j))/dy end doend do!------------------------------------------------------------------------------!内ú层?边?界?由?一?阶×迎?风?离?散Α?格?式?计?算?得?到?j=2do i=3,mx-1,1call upbound_v(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,vn)end doj=my-1do i=3,mx-1,1call upbound_v(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,vn)end doi=2do j=2,my-1,1call upbound_v(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,vn) end doi=mxdo j=2,my-1,1call upbound_v(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,vn) end do!----------------------------------------------------------------------!计?算?外猘层?边?界?do i=2,mx,1vn(i,1)=0.0vn(i,my)=0.0end dodo j=1,my,1vn(1,j)=-vn(2,j)vn(mx+1,j)=-vn(mx,j)end doreturnend subroutine!-----------------------------------------------------------------!定¨义?一?个?函ˉ数簓function alfa(x)implicit nonereal(kind=8) :: alfa,xif(x>0.0) thenalfa=1.0elsealfa=0.0end ifreturnend!利?用?一?阶×迎?风?格?式?求ó解a内ú层?边?界?上?的?速ù度è值μsubroutine upbound_u(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,un) implicit noneinteger :: i,j,mx,myreal(kind=8) :: dx,dy,dt,density,Viscosityreal(kind=8) :: u(mx,my+1),v(mx+1,my),p(mx+1,my+1)real(kind=8) :: un(mx,my+1)real(kind=8) :: df,dw,de,ds,dnreal(kind=8) :: aw,ae,as,an,apdw=Viscosity*dy/dxde=Viscosity*dy/dxds=Viscosity*dx/dydn=Viscosity*dx/dyaw=dw+max(0.5*density*(u(i-1,j)+u(i,j))*dy,0.0)ae=de+max(0.0,-0.5*density*(u(i,j)+u(i+1,j))*dy)as=ds+max(0.5*density*(v(i,j-1)+v(i+1,j-1))*dx,0.0)an=dn+max(0.0,-0.5*density*(v(i,j)+v(i+1,j))*dx)df=0.5*density*(u(i+1,j)-u(i-1,j))*dy+0.5*density*(v(i,j)+v(i+1,j)-v(i,j-1)-v(i+1,j-1)) *dxap=aw+ae+as+an+dfun(i,j)=u(i,j)+(dt/(dy*dx))*(-ap*u(i,j)+aw*u(i-1,j)+ae*u(i+1,j)+as*u(i,j-1)+an*u(i,j+1) )-dt*(p(i+1,j)-p(i,j))/dxreturnend subroutine!利?用?一?阶×迎?风?格?式?求ó解a内ú层?边?界?上?的?速ù度è值μsubroutine upbound_v(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,vn)implicit noneinteger :: i,j,mx,myreal(kind=8) :: dx,dy,dt,density,Viscosityreal(kind=8) :: u(mx,my+1),v(mx+1,my),p(mx+1,my+1)real(kind=8) :: vn(mx+1,my)real(kind=8) :: df,dw,de,ds,dnreal(kind=8) :: aw,ae,as,an,apdw=Viscosity*dy/dxde=Viscosity*dy/dxds=Viscosity*dx/dydn=Viscosity*dx/dyaw=dw+max(0.5*density*(u(i-1,j)+u(i-1,j+1))*dy,0.0)ae=de+max(0.0,-0.5*density*(u(i,j)+u(i,j+1))*dy)as=ds+max(0.5*density*(v(i,j-1)+v(i,j))*dx,0.0)an=dn+max(0.0,-0.5*density*(v(i,j)+v(i,j+1))*dx)df=0.5*density*(u(i,j)+u(i,j+1)-u(i-1,j)-u(i-1,j+1))*dy+0.5*density*(v(i,j+1)-v(i,j-1)) *dxap=aw+ae+as+an+dfvn(i,j)=v(i,j)+(dt/(dy*dx))*(-ap*v(i,j)+aw*v(i-1,j)+ae*v(i+1,j)+as*v(i,j-1)+an*v(i,j+1) )-dt*(p(i,j+1)-p(i,j))/dyreturnend subroutine!人?工¤压1缩?法ぁ?求ó解a下?一?时骸?刻ì的?压1强?值μsubroutine calp(un,vn,pn,mx,my,dx,dy,dt,c2,p)implicit noneinteger :: i,j,mx,myreal(kind=8) :: dx,dy,dt,c2real(kind=8) :: un(mx,my+1),vn(mx+1,my),pn(mx+1,my+1)real(kind=8) :: p(mx+1,my+1)!根ù据Yun和ívn求ó解a压1强?pn的?值μ!利?用?连?续?性?方?程ì求ó解a压1力do i=2,mx,1do j=2,my,1pn(i,j)=p(i,j)-dt*c2*((un(i,j)-un(i-1,j))/dx+(vn(i,j)-vn(i,j-1))/dy) end doend do!利?用?虚é拟a边?界?条?件t求ó解a压1力值μdo i=2,mx,1pn(i,1)=pn(i,2)pn(i,my+1)=pn(i,my)end dodo j=1,my+1,1pn(1,j)=pn(2,j)pn(mx+1,j)=pn(mx,j)end doreturnend subroutine。
计算流体大作业-基于C++的顶部驱动流模拟
-7-
计算流体力学大作业
第8页 共9页
图 6
Re=2000 时顶盖驱动流的流线图
图 7
Re=ห้องสมุดไป่ตู้000 时的顶盖驱动流流线图
-8-
计算流体力学大作业
第9页 共9页
图 8
Re=10000 时顶盖驱动流流线图
如图 6 至图 8 三幅图所示，基于 C/C++顶盖驱动流的流线图中，具有如下特 征： a) 方腔中心具有一个大漩涡； b) 方腔的左上角、左下角以及右下角具有次级漩涡； c) 随着雷诺数的增大，右下角的次级漩涡逐渐变小，左上角的次级漩涡逐 渐变大。
Re
ul Re ，其中顶盖流速 u = 0.1，l = 1
GAM 1 Re
4.程序组织
本文借助 D2Q9 格子波尔兹曼方法，在 Visual Studio 2012 平台上运用 C/C++语 言编程，运用有限差分方法，对上述传统问题进行了求解，并借助 Tecplot360 工具， 对所得数据进行后处理。所述程序组织如下图所示（完整程序参见所提交作业的文件 夹） ：
式中，p 为压力，v 为动力粘性系数。 根据上述分析，可以方便地设计出格子波尔兹曼算法。
（7）
3.2 顶部驱动流控制方程
方腔顶盖驱动流是考核程序的经典算例之一，其满足如下控制方程：
uu uv 2u 2u p 2 2 y x y x x uv vv 2v 2v p x 2 y 2 y y x
7.结论
本文借助 D2Q9 格子波尔兹曼方法， 在 Visual Studio 2012 平台上运用 C/C++ 语言编程，求解了顶部驱动流模型，运用 tecplot360 对所得数据进行处理，得 到了图 6 至图 8 三幅对应于不同雷诺数的流线图， 显示了该流动模型的典型现象， 验证了本文所述方法的正确性。
流场模拟方法
流场模拟方法流场模拟方法是一种重要的科学技术手段，用于研究和预测流体在各种条件下的运动和相互作用。
它在许多领域中都具有重要应用，如天气预报、风洞试验、环境工程和生物医学研究等。
流体力学是研究流体力学行为的学科，其中流场模拟方法是一个关键的研究领域。
流场模拟方法可以通过数学模型和计算机仿真来预测和分析流体流动的物理特性，从而为各种应用提供有效的解决方案。
流场模拟方法主要包括数值模拟和实验模拟两种。
数值模拟方法是通过建立数学模型和使用计算机算法来模拟流体运动。
这种方法的优点是可以准确预测流场的各种性质，如速度、压力、温度等，并能够在很短的时间内得到结果。
然而，数值模拟方法需要依赖复杂的数学模型和计算机算法，因此对计算资源要求高，而且模拟结果可能受到模型的假设和参数选择的影响。
实验模拟方法是通过设计和进行实验来模拟流体运动。
这种方法的优点是可以直接观测和测量流体的运动和相互作用，对结果的可信度高。
同时，实验模拟方法也能够提供丰富的数据来验证和改进数值模拟方法。
然而，实验模拟方法需要大量的设备和实验操作，并且受到实验条件和测量误差的限制。
在流场模拟方法中，数值模拟方法常用的技术包括有限元法、有限差分法和边界元法等。
这些技术通过对流体运动的偏微分方程进行离散化和求解，从而获得流场的数值解。
有限元法是一种广泛应用的数值模拟方法，它把流场划分为多个小单元，然后通过求解各单元上的方程来获得整个流场的数值解。
有限差分法是另一种常用的数值模拟方法，它将流场划分为网格点，在每个网格点上计算流体的变化量，然后通过迭代求解来获得整个流场的数值解。
边界元法是一种基于边界条件的数值模拟方法，它将流场划分为多个边界元，然后通过求解边界元上的方程来获得整个流场的数值解。
这些数值模拟方法都有各自的优点和适用范围，在具体应用中需要根据问题的复杂程度和计算资源的限制来选择合适的方法。
实验模拟方法中常用的技术包括风洞试验、流体力学实验和粒子图像测速法（PIV）等。
方腔顶盖驱动流数值模拟
方腔顶盖驱动流数值模拟张鑫（浙江理工大学 动力工程 2013G0502003）摘 要：在计算流体力学的研究中，通常要计算方腔驱动流问题来检验各种N-S 数值方法的有效性。
本文利用Fluent 软件对标准计算流体力学测试算例——方腔驱动流问题进行了模拟分析，其计算结果与文献中的标准解符合的比较好。
关键字：N-S 方程 方腔驱动流 Fluent 数值求解0引言流体流动的数值模拟广泛应用于气象、航天、机械、采矿等自然研究和工程计算的各个领域。
近年来，随着高性能计算与通信的迅速发展，针对流体流动的数值模拟以及求解相应Navier -Stokes 方程(简称N-S 方程)的高级算法研究现已成为目前国内外备受关注的热点和前沿课题。
Fluent 软件是用于模拟具有复杂外形的流体流动以及热传导的计算机程序，可以有效地模拟方腔驱动流问题，为计算流体力学的算法理论研究提供仿真参考。
高殿荣等学者采用液压冲击进行了分析；韩善玲等分析流体在空腔内的运动规律和物理机制，指出微小的凹凸是引起噪声的原因之一。
杨晶用Fluent 软件对方腔驱动流动进行了模拟分析，研究了不同雷诺数对计算结果的影响。
1模型介绍下图描述了本文所研究的物理模型，模型为边长等于0.1m 的正方形，上壁面为有一定速度的水，两侧壁面及地面均固定。
流体材料为水，密度为998.2kg/m3，黏度310005.1-⨯=u 。
2数值计算2.1、N-S 方程本文控制方程采用纳维司托克斯方程，纳维司托克斯方程是描述粘性不可压缩流体动量守恒的运动方程。
简称N-S 方程。
在直角坐标系中，可表达为如下所示： 连续方程：0=∂∂+∂∂yv x u 动量方程：）（yu x u x p y u v x u u 22221∂∂+∂∂+∂∂-=∂∂+∂∂υρ ）（yv x v x p y v v x v u 22221∂∂+∂∂+∂∂-=∂∂+∂∂υρ 2.2、网格划分及边界条件设置在gambit 软件中建立模型划分网络，由于模型几何形状比较规则，故全部采用四边形的的结构化网格，如下图所示。
EGKS 顶盖驱动方腔流动数值模拟研究
EGKS 顶盖驱动方腔流动数值模拟研究作者：杨约翰
来源：《科学与财富》2020年第20期
摘要：稀薄气体条件下，气体分子的宏观性质会受到气体分子的间断粒子效应的影响，连续介质模型不能完全反映物理实际。
建立在连续介质模型基础上的流体力學控制方程 NS 方程就不能适用于此。
本文选用 EGKS 模型的 Blotzmann 方程作为其控制方程进行了数值模拟，对不同 Kn 数下顶盖驱动方腔流动的流动特性进行了研究。
关键词：稀薄气体;NS 方程;Blotzmann 方程;EGKS 模型;顶盖驱动方腔流动。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
方腔顶盖驱动流流场数值预测摘要：本文分别采用一阶迎风格式（FUD）、中心差分格式(CD)和乘方格式(PLD)计算方腔顶盖驱动流，计算结果同Ghia et al结果进行比较。
由计算结果可得出，一阶迎风引起的假扩散最大，计算结果偏离基准解最远，中心差分格式和乘方格式同基准解已经非常接近。
但中心差分格式不稳定，不易收敛。
网格数变化也会对结果产生影响，网格划分越多，计算结果与基准解越接近，而计算的时效性越差，所以在划分网格时，我们需要综合考虑其准确性和时效性，选用合理网格数。
关键字：一阶迎风格式，中心差分格式，乘方格式，网格数The prediction of flow field in the flow in driven cavity Abstract:In this paper, the three discrete formats of the equation convection (PLD, FUD and CD) was used to calculation the flow field in the flow in driven cavity. Through the compared with Ghia et al, we found that the false diffusion is the largest caused by the FUD, and the deviation of the calculation results from the exact solution, CD is the least , PLD come next and FUD is the largest. But CD is instability, it’s difficult converg ence. The changes of grid number will have an impact on the results. By the analysis, the more grid, the closer of the calculated results with the exact solution, and the worse of the calculated timeliness, so meshing, we need consideration of it’s accurac y and timeliness, to get a reasonable number of grid.Key words: FUD ，CD，PLD, the number of grid引言对流-扩散方程离散格式的稳定性与准确性一直是数值传热学中的一个重要问题，而对流-扩散方程的离散关键在于对流项的离散。
对流项常见的离散格式有乘方格式（PLD），一阶迎风格式（FUD），中心差分格式（CD）,这三种格式在计算精度和计算时效上各有优缺点。
方腔顶盖驱动流是考核程序的经典算例之一，本文就以上三种格式在雷诺数分别为100、400、1000、3200的情况下对方腔顶盖驱动流流场进行数值预测，并将其计算结果与Ghia et al结果进行对比分析。
1. 控制方程u=0v=0llu=1,v=0u=0v=0u=0,v=0xy()()()()22222222uu uv p u u x y x x y uv vv p v v xy y x y ρηρη∂∂⎡⎤⎛⎫∂∂∂+=-++ ⎪⎢⎥∂∂∂∂∂⎝⎭⎣⎦∂∂⎡⎤⎛⎫∂∂∂+=-++ ⎪⎢⎥∂∂∂∂∂⎝⎭⎣⎦()()()()2222222211uu uv puu x y x x y uv vv pvv xyy x y ηρρηρρ∂∂⎛⎫∂∂∂+=-++ ⎪∂∂∂∂∂⎝⎭∂∂⎛⎫∂∂∂+=-++ ⎪∂∂∂∂∂⎝⎭Re Re ul ρρηη=⇒=，其中顶盖流速u = 1，l = 1 图1 方腔顶盖驱动流1G A M R e=2.对流项离散格式所谓对流项离散格式，就是指控制容积界面上被求物理量的插值方式[2]， 恒取上游节点的值，即：1,0E e ea P D ∆=+⎡-⎤⎣⎦而中心差分则取上、下游节点的算术平均值，即：112E eea P D ∆=-对于乘方格式，则有：()50,10.10,Eee ea P P D ∆∆⎡⎤=-+⎡-⎤⎣⎦⎢⎥⎣⎦3. 格式对比分析3.1 三种格式计算结果与Ghia et al 结果对比分析图2给出了在129⨯129均匀网格下采用三种格式的计算结果[3]与Ghia et al 结果进行的对比分析。
从图中可以看出，一阶迎风(FUD)引起的假扩散最大，计算结果偏离基准解最远，且Re 数越大，偏差越大；而乘方格式（PLD ）、中心差分格式(CD)的计算结果都已经逐渐靠近基准解。
值得一提的是，对于中心差分格式(CD)在Re=3200，松弛因子为0.5时计算结果是发散的，图中给出的是松弛因子降为0.2时的计算结果。
并由图中可知，中心差分格式(CD)与乘方格式（PLD ）0.00.20.40.60.81.0-0.4-0.20.00.20.40.60.81.0 1.2U (m /s )y Re=100Ghia FUD CD PLD0.00.20.40.60.8 1.0-0.4-0.20.00.20.40.60.81.0 1.2Ghia FUD CD PLDU (m /s )yRe=4000.00.20.40.60.81.0-0.4-0.20.00.20.40.60.81.0 1.2Ghia FUD CD PLDU (m /s )yRe=10000.00.20.40.60.81.0-0.6-0.4-0.20.00.20.40.60.81.01.2GhiaFUD CD PLDU /(m/s )yRe=32000.00.20.40.60.8 1.0-0.3-0.2-0.10.00.10.2Ghia FUD CD PLDV (m /s )xRe=1000.00.20.40.60.8 1.0-0.5-0.4-0.3-0.2-0.10.00.10.20.30.4Ghia FUD CD PLDV (m /s )xRe=4000.00.20.40.60.8 1.0-0.6-0.4-0.20.00.20.4Ghia FUD CD PLDV (m /s )xRe=10000.00.20.40.60.8 1.0-0.6-0.4-0.20.00.20.40.6Ghia FUD CD PLDV (m /s )xRe=3200图2三种格式与Ghia 计算结果比较相比在Re=3200时更加接近基准解，由此可知若中心差分格式(CD)可以取的收敛的解，则其精度更高。
而图中其他格式松弛因子均为0.5，同样中心差分格式(CD)在雷诺数为100、400和1000时松弛因子也为0.5. 由此可得结论中心差分格式(CD)具有较高的准确性但其不具有稳定性。
3.2 三种格式流场图和压力场的比较分析图3、4给出了在129 129均匀网格下，Re=3200时采用三种格式的计算结果所绘制的流场图和压力场。
从图中可以看出，乘方格式（PLD）和中心差分格式(CD)的流场图较为相似，一阶迎风(FUD)的流场与前两者相比有较大差异，由前述可知，此为一阶迎风(FUD)引起的假扩散。
另外，图中所示为中心差分格式(CD)松弛因子降为0.2时的计算结果，因其精度高于乘方格式（PLD），从图中可以看出其流场模拟具有更高的准确性。
PLD-0.5 CD-0.2 FUD-0.5图3三种格式流场图PLD-0.5 CD-0.2 FUD-0.5图4三种格式压力场图3.3 三种格式的计算时效性从表1计算时间来看，随着雷诺数和网格数的增大，计算时间明显增多，而三种格式之间并无明显规律。
由表2可得，中心差分格式(CD)具有不稳定性，随着松弛因子的减小，最终可以获得收敛的解，而计算时间和迭代次数也明显增多。
表1 方腔顶盖驱动流（relax=0.5，SMAX=1.0E-8）Re网格数 时间s迭代次数 FUD CD PLD FUD CD PLD 3200 414 发散 4 967 发散 1026 81 43 发散 46 2604 发散 2827 129174 发散 263 4627 发散 4852 1000 129 163 113 154 3303 2820 2946 400 158 162 154 3210 3155 3135 10080104123196620852086表2 中心差分格式(CD)网格数 Re=3200松弛因子 0.5 0.2 0.02 41 迭代次数 发散 发散 7522 时间s 33 81 迭代次数 发散 5675 19009 时间s 98 306 129迭代次数 发散7820 7505 时间s2872763.4 网格数变化对计算结果的影响图5为不同网格数下计算结果与基准解进行比较分析，由图中可知，网格划分越多，计算结果与基准解越接近，而网格划分的越密，计算的时效性越差，迭代的次数也相应变多，所以在划分网格时，我们需要综合考虑其准确性和时效性。
0.00.20.40.60.81.0-0.6-0.4-0.20.00.20.40.60.81.01.2Ghia 41x41 81x81 129x129U (m /s )yRe=3200PLD 0.00.20.40.60.8 1.0-0.6-0.4-0.20.00.20.4Ghia41x41 81x81 129x129V (m/s )xRe=1000FUD图5网格数变化对计算结果的影响4. 结论本文以一阶迎风格式（FUD），中心差分格式（CD）和乘方格式（PLD），分别对方腔顶盖驱动流的流场进行数值预测，并将计算结果同Ghia et al结果进行比较，考察了三种格式的计算精度与计算时间。
计算表明，一阶迎风引起的假扩散最大，计算结果偏离基准解最远，中心差分格式和乘方格式同基准解已经非常接近，并且中心差分格式(CD)与乘方格式（PLD）相比在Re=3200时更加接近基准解，其精度更高。
但中心差分格式的不稳定性，若要获得收敛的解就必须减小松弛因子，随之而来的便是更多的计算时间和迭代次数。
综合考虑各种因素后：这三种格式的计算精度和计算时间都不甚合理，有待于使用更高精度和时效性的格式对该问题进行计算。
网格数变化也会对计算结果产生影响，网格划分越多，计算结果与基准解越接近，而网格划分的越密，时效性越差，迭代次数也变多，所以在划分网格时，我们需要综合准确性和时效性，定出合理的网格数。
5. 参考文献[1] Ghia U, Ghia K N and Shin C T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multi-grid method[J]. J Comput Phys, 1982. 48: 387-411[2] 陶文铨编著. 数值传热学（第二版），西安：西安交通大学出版社，2001[3] 传热与流体流动的数值计算课程教学程序，EXAMPLES。