方腔顶盖驱动流流场数值预测
方腔顶盖驱动流动
一、二、问题描述方腔顶盖驱动流动如图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 模型;顶盖驱动方腔流动。
N-S方程基于投影法的特征线算子分裂有限元求解
N-S方程基于投影法的特征线算子分裂有限元求解水庆象;王大国【摘要】提出了一种求解非定常不可压缩纳维-斯托克斯方程(N-S方程)的新型有限元法:基于投影法的特征线算子分裂有限元法.在每一个时间层上将N-S方程分裂成扩散项、对流项、压力修正项.对流项采用多步显式格式,且在每一个对流子时间步内采用更加精确的显式特征线-伽辽金法进行时间离散,空间离散采用标准伽辽金法.应用此算法对平面泊肃叶流、方腔流和圆柱绕流进行数值模拟,所得结果与基准解符合良好.尤其对于Re=10000的方腔流,给出了方腔中分离涡发展和运动的计算结果,并发现在该雷诺数下存在周期解,表明该算法能较好地模拟流体流动中的小尺度物理量以及流场中分离涡的运动.【期刊名称】《力学学报》【年(卷),期】2014(046)003【总页数】13页(P369-381)【关键词】非定常不可压纳维-斯托克斯方程;投影法;特征线算子分裂有限元法;多步格式【作者】水庆象;王大国【作者单位】西南科技大学环境与资源学院,绵阳621010;西南科技大学环境与资源学院,绵阳621010【正文语种】中文【中图分类】O357.1有限元法被广泛应用于流体动力学问题的求解，并表现出良好的几何边界和边界条件适用性.而将标准伽辽金有限元法用于非定常不可压纳维--斯托克斯方程的求解时，会出现数值振荡.主要原因有两点[1]：一是由于压力变量不出现在连续方程中，离散过程中速度场和压力场有限元插值函数的组合不恰当，使得Ladyˇzenskaja--Babuˇska--Brezzi(LBB)条件不能满足，从而引起压力场的数值振荡；二是动量方程中存在非线性的对流项，当雷诺数较高对流占优时，标准伽辽金有限元求解将导致数值振荡.为避免速度--压力插值函数不满足LBB条件而导致压力场的数值振荡，Chorin[2]提出了投影法.投影法的原理是通过引入一个无需满足散度为零的中间变量对压力和速度解耦求解，速度和压力解耦后无需要求速度--压力的插值函数满足LBB条件，并且能够高效地计算非定常问题[3].Kim和Moin[4]提出了具有二阶精度的投影法；Donea等[5]将投影法从有限差分领域推广至有限元领域；Brown等[6]对于不同的投影法作了具体分析并提出了改进方法；刘淼儿等[7]在对投影法研究进展进行综述的基础上，发展了在时间上具有三阶精度的投影法.针对对流占优导致的数值解振荡，学者们发展了多种稳定化有限元法解决该问题.目前应用较为广泛的有流线迎风彼得洛夫--伽辽金(streamline upwind Petrov--Galerkin,SUPG 法)[8]、泰勒 --伽辽金法[9]、特征线--伽辽金法[10]等.其中，特征线--伽辽金法沿特征线对时间进行离散，以理性的形式引入稳定项，具有明确的物理意义，可以有效克服数值振荡，从而保证数值解的稳定，空间上采用标准伽辽金法进行离散.Wang等[11]将泰勒展开引入到特征线--伽辽金法并结合算子分裂法的优点提出基于特征线的算子分裂(characteristic-based operator splitting algorithm,CBOS)有限元法.该方法将纳维--斯托克斯方程分成扩散项和对流项，对流项时间离散借鉴了基于特征线分裂(characterictic-based split，CBS)算法[12]，沿特征线进行离散给出合理的平衡耗散项，从而避免了彼得洛夫--伽辽金法等其他有限元法修正权函数的困难.但是该算法在处理扩散项时对速度、压力进行耦合求解，为了使速度--压力的插值函数满足LBB条件，一般要求压力插值函数比速度插值函数低一阶，这使得研究者不能任意选取速度--压力插值函数；对流项显式求解，采用显式格式对时间步长的要求较为苛刻，计算成本较高；同时为了获得完全显式格式采用的近似关系在时间步长稍大时会带来一定的误差.本文在文献[2,11]的基础上，提出一种求解不可压黏性流问题的新型有限元法：基于投影法的特征线算子分裂有限元法，并实现了该算法有限元数值计算程序的开发.为验证算法的有效性，对平面泊肃叶流、方腔流和圆柱绕流进行数值模拟，并将计算结果与已有结果进行了对比分析.二维非定常不可压缩流体的N--S方程无量纲形式表示为[13]式中，Ω和[0,T]分别为空间域和时间域，i,j=1,2，)，u为水平方向速度，v为垂直方向速度，p为压力，t为时间，为雷诺数(其中，µ为黏性系数，ρ为流体密度，U为特征速度，L为特征长度)，f1为水平方向外力，f2为垂直方向外力，(x1,x2)=(x,y).2.1 N--S方程的分裂基于投影法和CBOS有限元法[11]的思想，把纳维--斯托克斯方程分裂成扩散项对流项压力修正项对式(4)两边取散度，结合式(5)可得压力泊松方程式(2)～式(6)中，和pn分别为n时刻的速度场和压力场的值，为扩散项(2)在n+1时刻的解，同时也为对流项(3)在n+1时刻的初值，为对流项(3)在n+1时刻的解，同时也为压力修正项(4)的初值，和和pn+1分别为n+1时刻的速度场和压力场的值，Δt为整体时间步长.2.2 对流项沿特征线的显式时间离散对流项(3)为双曲型方程，其特征线方程为式中，xj为特征线上位置的分量，uj为特征速度分量.式(7)用差分形式可写为式中，为在时间间隔Δt内特征速度分量uj的平均值.沿着特征线，式(3)可写为比较式(3)与式(9),沿着特征线方向，对流项的非线性部分消失，空间离散可采用标准伽辽金有限元法，从而避免SUPG和GLS[1,12]等有限元法合理选择权函数的困难.式(9)在时间域的差分形式为将式(8)代入式(11)，得式中，k=1,2.为了得到完全显式的时间离散格式，一般取，当时间步较大时会带来一定误差.为了获得更为精确的完全显式离散格式，对在处进行泰勒展开[14]，得将式(14)代入式(13)，并略去O(Δt3)，有与CBOS有限元法[11]相比，式(15)中的右边第2项即为考虑时间高阶量影响的附加稳定项.2.3 对流项的多步格式显式格式是条件稳定的，受到了时间步长的严重束缚.本文采用多步格式[15]处理对流项以降低对整体时间步长的要求，提高算法的稳定性，即在每一个整体时间步内将对流项分为h步进行求解，则求解对流项的时间步长.根据式(15)，对流项在第l子时间步内的求解公式为将式(17)代入式(16)，有3.1 扩散项有限元离散对扩散项(2)采用向后差分格式进行时间离散由伽辽金法[16]建立式(19)的弱形式，并对黏性项分部积分,可得式中，Γ为区域边界，δui为虚位移.按有限元求解方法，将计算域划分为E个单元并将一典型单元记为Ω(e)，取典型单元分析式中，Γ(e)为典型单元边界，取典型单元为四节点四边形单元，令式中，α=1,2,···,m，m=4为插值单元内速度节点个数.将式(22)代人式(21)，并令，得单元有限元方程式中，β=1,2,···,m,式(24)中，δij为置换算子，δijuiβ=ujβ.在整个计算域合成可得总体有限元方程s,q为速度总体节点号.解式(28)可得速度uin+θ1.3.2 对流项有限元离散同扩散项处理方法一样，由标准伽辽金法建立式 (18)的弱形式，并简记上标，，得对上式最后一项分部积分，忽略边界项的影响，得由式(22)可得式(31)中，Γ,η=1,2,3,4.按有限元求解方法得单元有限元方程式中在整个计算域合成可得总体有限元方程解式(34)可得速度ug(l+1)i.3.3 压力泊松方程有限元离散由标准伽辽金法建立压力泊松方程(6)的弱形式式中，δp为虚位移.取典型单元为四节点四边形单元，令得单元有限元方程式中在整个计算域合成可得总体有限元方程解式(42)可得压力pn+1.3.4速度校正项有限元离散式(4)整理可得由标准伽辽金法建立方程(43)的弱形式依据有限元求解方法得单元有限元方程在整个计算域合成可得总体有限元方程解式(46)可得速度un+1i.当n时刻的值已知，求n+1时刻的值时，纳维--斯托克斯方程的求解过程为(1)以n时刻的速度场和压力场 pn作为初值，根据式(2)求解速度场；(2)以为初值，由式(3)求解速度场.在求解式(3)时，采用了多步显示格式求解，具体求解流程如图1所示；(3)已知，求解压力泊松方程式(6)得压力场pn+1；(4)根据式(4)得到修正后的速度，便完成了由n时间步到n+1时间步的求解；(5)转到下一时刻，重复步骤(1)～(4).5.1 平面泊肃叶流平面泊肃叶流是少数存在理论解的流动问题之一，其几何模型及边界条件如图2所示.计算模型中y方向取特征长度L=1m，腔体长高比取6:1，特征速度U=1m/s，液体密度ρ=1.0kg/m3.入口处的无因次速度按抛物线分布为u=6y(1−y),v=0，固壁面上满足不可滑移条件u=0,v=0，出口处压力p=0.整个计算区域剖分为10×60个四节点四边形单元.平面泊肃叶流的压力理论解为[17]取黏性系数µ=0.01Pa·s，则Re=100.在计算过程中，判断流场稳定的标准为检验相邻两个时间步变量的变化量[18]式中，φ为速度、压力的任意量，本文选取φ=u，i为从1变化到整个计算域内节点总数.为比较对流项多步格式采用不同步数h的不同数值表现，采用本文算法基于不同的时间步长Δt和不同步数h对Re=100的平面泊肃叶流进行数值模拟.水平中轴线(x,0.5m)上压力的计算结果与理论解的相对误差定义为[19]式中，p∗为本文计算结果，p∗∗为理论解，Num为中轴线(x,0.5m)上节点总数.表1列出了相对误差η随不同时间步长Δt和不同h的变化.由表1可知，在Δt=0.02情况下，取h=1和h=2时，本文算法不能收敛；在Δt=0.01情况下，取h=1时，本文算法不能收敛；在Δt=0.005情况下，无论h如何取值本文算法均收敛，并且相对误差η随着h的增加而减小.由以上分析可知对流项采用多步格式既能降低对整体时间步长的要求，提高算法稳定性，又可提高算法精度.5.2 方腔流顶盖驱动方腔流常被用来检验求解黏性不可压缩纳维--斯托克斯方程数值算法的可靠性，其几何模型及边界条件如图3所示.取特征长度L=1m，特征速度U=1m/s，液体密度ρ=1.0kg/m3，通过调整黏性系数得到不同的雷诺数.因此，方腔无因次的边长为1，顶部施加无因次水平速度为1，垂向速度设为0，其他3个边界均为不可滑移固壁边界条件u=0,v=0，左下角指定相对压力p=0，对边界附近的网格进行加密.本节对方腔流进行数值计算时取h=10.5.2.1 稳定层流解的分析Liu等[20]和Lin等[21]指出当Re≤5000时顶盖驱动方腔流内部存在稳定的层流解，采用本文算法在60×60四节点四边形网格下对Re=100,400,1000,3200,5000的方腔流进行数值模拟.在计算过程中，判断流场稳定的标准如式(48).图4给出了不同雷诺数下水平速度u沿垂直中线、垂向速度v沿水平中线分布结果.图中实线表示本文计算结果；点表示Ghia等[22]的计算结果，其结果被视为标准解.由图4可见本文计算结果与Ghia等计算结果十分吻合.在Ghia等模型中，当Re=100, 400,1000,3200时，网格数为129×129；Re=5000时,网格数为257×257，而本文采用网格数仅为60×60，远少于Ghia等模型中采用的单元数.以上分析表明本文算法在采用较少网格的情况下也能达到较高的精度.图5给出了Re从100到5000各状态下压力等值线图(上)和涡量图(下)，方腔流中的压力等值线较为光滑，没有出现振荡.图6给出了Re从100到5000各状态下流场稳定后流线图，并包括一些局部流场的放大图.从图6(a)中可以看到Re=100时方腔底部两角各产生一个二次涡，右边的二次涡较大.由图6(b)可见，随着Re增加方腔底部两个二次涡也随着增大.当Re 增至1000时，方腔左下角的二次涡有明显的增大但仍保持一个二次涡，如图6(c)所示.从图6(d)中可以看到Re=3200时，左下角二次涡进一步增大，右下角新产生一个次级二次分离涡，并且在方腔左固壁的上半部也分离出一个二次涡.当Re=50000时，此时流场与Re=3200时的结果相似，只是所有分离涡大些，尤其是右下角产生的次级二次分离涡变得很明显，如图6(e)所示.这些图像与Ghia等[22]计算所得流动图像基本相同，说明本文算法能够很好地模拟流动中的小尺度物理量.5.2.2 周期解的分析采用本文算法对Re=10000方腔流进行数值模拟，网格数为100×100.图7给出了Re=10000在点(0.125,0.9375)处水平方向速度u(左图)和垂直方向速度v(右图)随时间的变化规律.由图可知，Re=10000时速度随时间作周期性的变化，对图7中的数据进行快速傅里叶变换可得周期变化的频率为0.6254，Bruneau和Sadd[23]在网格数为512×512时的周期频率为0.61.图8给出了Re=10000方腔流流场中流线随时间变化规律，每两幅图之间的时间间隔为0.25.比较图8(a)和图8(b)可知随时间发展左下角的二次涡向左收缩且在其右边开始产生一个新二次涡；同时，左上角的二次涡出现了向上收缩的现象.在图8(c)中，左下角新产生的二次涡更为明显，而原二次涡则继续向左收缩；左上角的二次涡继续向上收缩，其下部固壁上产生了一个新二次涡.从图8(d)中可知左下角新产生的二次涡进一步增大，而原二次涡继续向左收缩且脱离方腔底壁；同时，左上角的二次涡下部固壁上新产生的二次涡进一步增大.在图8(e)中，左下角新产生的二次涡比原来的二次涡略大，左上角二次涡进一步向上收缩，其下部新产生的二次涡更为明显.在图8(f)中，左下角新产生的二次涡相对原二次涡大进一步增大，同时左上角新产生的二次涡开始与原二次涡合并.在图8(g)中左下角与左上角新旧二次涡都已合并成一个二次涡.图8(a)与图8(g)两幅图像之间的时间间隔为1.5，在文献[20]中指出Re=10000的方腔流流动图像变化的周期为1.59，所以图8(a)与图8(g)较为接近.Bruneau和Sadd[23]采用512×512的网格也得到了相同的流线周期性变化规律.以上分析结果表明，本文算法能够在较少网格下很好地模拟高Re数流场中分离涡的运动情况.5.2.3 求解顺序的影响分析基于Re=1000的方腔流从收敛速度和计算精度两个方面比较了算法1和算法2，其中算法1为先求解扩散项后求解对流项，算法2为先求解对流项后求解扩散项. 图9给出了两种算法计算过程中的误差ε与计算步数的对应关系.图中实线为算法1模拟结果；虚线为算法2模拟结果.由图9可见算法1与算法2在收敛速度上是基本一致的.图10给出了两种算法计算所得水平速度u沿垂直中线、垂向速度v沿水平中线分布结果.图中实线表示算法1模拟结果；虚线表示算法2模拟结果；点表示Ghia等[22]的计算结果.由图10可见算法1的精度优于算法2，即本文所采用的先求解扩散项后求解对流项的策略在计算精度上占优.5.3 圆柱绕流计算区域尺度在主流方向上取为30D，其中圆柱上游分配10D，横向尺寸为18D，圆柱直径D= 1m为特征长度.整个计算域划分为11860个四节点四边形单元，共有12114个节点，其中圆柱表面分布128个网格.入口处指定沿水平方向的均匀来流U=1m/s，垂向速度0；侧壁采用可滑移边界条件；圆柱表面为不可滑移边界条件；出口处为自由出流边界，并指定压力为p=0.液体密度ρ=1.0kg/m3，通过调整黏性系数得到不同的雷诺数.本节对圆柱绕流进行数值计算时取h=10.图11给出了Re=200时阻力系数和升力系数的时程曲线.表2列出了Re=100,200时本文计算的阻力系数Cd，升力系数Cl和斯特劳哈数St与相应文献[24-26]结果的比较.由表2可见，本文结果与其他结果接近，进一步说明了本文算法是准确和可靠的.在 CPU为 i7-2.93GHz且 RAM为 3.49G的PC机上，分别采用本文算法和Wang等[11]提出的CBOS有限元法对Re=5000方腔流进行数值模拟.6.1 计算效率的比较投影法直接对纳维--斯托克斯方程中的速度和压力进行解耦，解耦后得到几个简单方程分别对其进行求解，从而得到速度和压力值.本文算法基于投影法将速度和压力进行解耦，而Wang等提出的CBOS有限元法在处理扩散项时仍是将速度和压力耦合求解，因此本文算法具有更高计算效率.表3列出两种算法的计算效率，由表3可知本文算法计算速度是CBOS有限元法[11]计算速度的1.6倍.6.2 计算精度的比较本文计算结果与Ghia等[22]计算结果的平均相对误差绝对值，uij为本文计算结果，为标准解，N′为Ghia等所给出的速度值的个数，N′=17)为3.218%，采用CBOS 有限元法计算得到的值为5.2%.因此，本文算法与CBOS有限元法[11]相比具有更高的计算精度，其原因主要是对流项采用了多步显式格式.此外，本文算法允许压力和速度采用任意阶次的同阶插值函数，更易于程序的实现.(1)本文提出了一种求解非定常不可压纳维--斯托克斯方程新型有限元方法：基于投影法的特征线算子分裂有限元法.该方法将投影法和基于特征线的算子分裂有限元法相结合，在每一个时间层上将纳维--斯托克斯方程分裂成扩散项、对流项、压力修正项.扩散项时间离散采用向后差分格式，空间离散采用标准伽辽金有限元法，隐式求解；对流项采用多步显式格式，且在每一个对流子时间步内采用更加精确的显式特征线--伽辽金法进行时间离散，空间离散采用标准伽辽金法，通过平面泊肃叶流模拟结果表明对流项采用多步格式既降低了对整体时间步长的要求，提高算法稳定性，又能提高算法精度；压力通过求解压力泊松方程获得，压力泊松方程结合连续方程推导得到，对压力泊松方程用标准伽辽金法进行空间离散；得到压力值后对速度进行修正.速度和压力解耦后无需要求速度--压力的插值函数满足LBB条件，大大提高求解效率，更易于程序的实现.(2)对不同雷诺数下的方腔流进行数值模拟，模拟结果与标准解符合很好，表明该算法具有较高的精度和稳定性；相对文献[11]中的算法，本文算法在计算精度和计算效率上也有较大的提高.进一步模拟了Re=10000的方腔流，在该雷诺数下存在非定常周期解，并给出方腔中分离涡的周期运动，表明该算法能较好地模拟流体流动中的小尺度物理量，可以将该算法用于研究在高雷诺数下流动从层流向湍流转变过程中分离涡运动的规律.(3)对于Re=100和200的圆柱绕流，计算所得的结果如阻力系数、升力系数、斯特劳哈数等与已有数据也较为接近.以上的对比结果表明本文建立的算法用于模拟圆柱层流绕流是准确的和可靠的，为进一步求解复杂的非定常流动问题奠定了基础.1 Hannani SK,Stanislas M,Dupont P.Incompressible Navier--Stokes equations with SUPG and GLS formulations—A comparisonputer Methods in Applied Mechanics and Engineering,1995, 124:153-1702 Chorin AJ.Numerical solution of the Navier--Stokesequations.Mathematics and Computation,1968,22:745-7623王坪,田振夫.基于投影法求解不可压缩流的高精度紧致格式.工程数学学报,2010,27(1):25-232(Wang Ping,Tian Zhenfu. A projection-based high-order compact scheme for solving incompressible fl w.Chinese Journal of Engineering Mathematics,2010, 27(1):25-232(in Chinese))4 Kim J,Moin P.Application of a fractional-step method to incompressible Navier--Stokes equations.Journal of Computational Physics,1985,59:308-3235 Donea J,Giulianli S,Laval H.Finite element solution of the unsteady Navier--Stokes equations by a fractional step puter Methods in Applied Mechanics and Engineering,1982,30(1): 53-736 Brown DL,Cortze R,Minion ML.Accurate projection methods for incompressible Navier--Stokes equations.Journal of Computational Physics,2001,168:464-4997刘淼儿,任玉新,张涵信.求解不可压Navier--Stokes方程的三阶精度投影方法.清华大学学报,2005,45(2):285-288(Liu Miaoer, Ren Yuxin,Zhang Hanxin. Third-order projection method for solving the incompressible Navier--Stokes equations.Journal of Tsinghua University(Science and Technology),2005,45(2):285-288(in Chinese))8 Ganesan S.An operator-splitting Galerkin/SUPG finit elment methodforpopulationbalanceequations:Stabilityandconvergence.ESAIM:Ma thematical Modelling and Numerical Analysis,2012, 46(6):1447-14659 Wang X,Ouyang J.Time-related element-free Taylor--Galerkin method with non-splitting decoupling process for incompressible steady flw.International Journal for Numerical Method in Fluids, 2012,68(7):839-855 10 Ding DY,Wu SQ.Direct numerical simulation of turbulent fl w over backward-facing at high Reynolds numerical.Science China Technological Science,2012,55(11):3213-322211 Wang DG,Wang HJ,Xiong JH,et al.Characteristic-based operatorsplittingfinit elementmethodforNavier--Stokesequations.Science China Technological Science,2011,54(8):2157-216612 Malan AG,Lewis RW.An artificia compressibility CBS method for modellingheattransferandflui fl winheterogeneousporousmaterials.International Journal for Numerical Methods in Engineering, 2011,87(1-5):412-42313林建忠,阮晓东,陈邦国等.流体力学.北京:清华大学出版社，2005(Lin Jianzhong,Ruan Xiaodong,Chen Bangguo,et al.FluidMechanics.Beijing:Tsinghua University Press,2005(in Chinese))14韩向科,钱若军,袁行飞等.改进的基于特征线理论的流体力学有限元法.西安交通大学学报,2011,45(7):112-117(Han Xiangke, Qian Ruojun,Yuan Xingfei,et al.Improved characteristic-based split finit element method for flui dynamics.Journal of Xian Jiaotong University,2011,45(7):112-117(in Chinese))15 Hundsdorfer W,Mozartova A,Spijker MN.Stepsize restrictions for boundedness and monotonicity of multistep methods.Journal of Scientifi Computing,2012,50(2):265-28616章本照.流体力学中的有限元方法.北京:机械工业出版社,1984 (Zhang Benzhao,The Finite Element Method in Fluid Dynamics. Beijing:ChinaMachine Press,1984(in Chinese))17 Li X,Han X.An iterative stabilized fractional step algorithm for numerical solution of incompressible N--S equations.International Journal for Numerical Methods in Fluids,2005,49(4):395-41618 Nithiarasu P.An efficient artificia compressibility(AC)scheme based on the characteristic based split(CBS)method for incompressible flws.International Journal for Numerical Methods inEngineering,2003,66(10):1815-184519 Li X,Duan Q.Meshfree iterative stabilized Taylor--Galerkin and characteristic-based split(CBS)algorithms for incompressible N--S puter Methods in Applied Mechanics and Engineering,2006,195(44):6125-614520 Liu H,Fu DX,Ma YW.Upwind compact method and dirct numerical simulation of driven fl w in a square cavity.Science in China(SeriesA),1993,23(6):657-66521 Lin LS,Chang HW,Lin CA.Multi relaxation time lattice Boltzmann simulations of transition in deep 2D lid driven cavity usingputers&Fluids,2013,80:381-38722 Ghia U,Ghia KN,Shin CT.High-Resolutions for incompressible fl w using the Navier--Stokes equations and a multigrid method.Journal of Computational Physics,1982,48(3):387-41123 Bruneau CH,Saad M.The 2D lid-driven cavity problemputer and Fluids,2006,35(3):326-34824詹昊,李万平,方秦汉等.不同雷诺数下圆柱绕流仿真计算.武汉理工大学学报,2008,30(12):129-132(Zhan Hao,Li Wanping,Fang Qinhan,et al.Numerical simulation of the fl w around a circularcylinderatvariesReynoldsnumber.JournalofWuhanUniversity of Technology,2008,30(12):129-132(in Chinese))25 Xu S,Wang ZJ.An immersed interface method for simulating the interaction of a flui with moving boundaries.Journal of Computa-tional Physics,2006,216(2):454-49326 Harichandan AB,Roy A.Numerical investigation of low Reynolds number fl w past two and three circular cylinder using unstructured grid CFR scheme.International Journal of Heat and Fluid Flow, 2010,31:154-171。
顶盖驱动流 谱方法
顶盖驱动流谱方法
顶盖驱动流谱方法（Top Hat driven flow spectroscopy method）是一种常用于表征复杂流体动力学行为的流体力学方法。
顶盖驱动流谱方法通过在流体上方施加一个顶盖运动，并通过监测流体中的扩散和剪切变形来评估流体的流动性质。
顶盖驱动流谱方法的基本原理是利用顶盖的运动来引起流体的流动。
通过在流体上方施加一个顶盖的运动，可以在流体中引入额外的扩散和剪切变形。
顶盖的运动可以通过不同的方式实现，如旋转、振荡或往复移动。
同时，在顶盖运动的影响下，流体中的粒子或标记物会发生位移和旋转，这些运动可以通过光学、电学或其他检测方法进行监测和分析。
顶盖驱动流谱方法可以用于研究不同流体性质的变化，如黏度、流变特性和颗粒浓度等。
通过测量扩散和剪切变形的程度，可以得到关于流体动力学行为的定量信息。
这种方法在生命科学领域的应用也很广泛，可以用于研究细胞运动、颗粒流和微流体中的流动等。
与传统的流体力学方法相比，顶盖驱动流谱方法具有一些优势。
首先，它可以对复杂、非牛顿流体的性质进行研究，而传统方法往往只适用于牛顿流体。
其次，顶盖驱动流谱方法可以用于研究微观尺度下的流体行为，可以探索纳米和微米尺度上的流动现象。
此外，这种方法还可以与其他技术相结合，如荧光标记、微流控等，进一步拓展其应用范围。
总之，顶盖驱动流谱方法是一种有效的流体力学方法，可以用
于研究复杂流体的动力学行为。
其原理简单，应用范围广泛，有望在生物医学和纳米技术等领域发挥重要作用。
基于数值模拟的混流式水轮机顶盖内流场分析
基于数值模拟的混流式水轮机顶盖内流场分析发表时间：2018-04-20T10:32:08.350Z 来源：《建筑学研究前沿》2017年第32期作者：李万[导读] 混流式水轮机的应用最为广泛，这种机组形式具有结构紧凑，运行可靠。
哈尔滨电机厂有限责任公司黑龙江哈尔滨 150040 摘要：本文以某高水头电站为研究对象，采用数值模拟方法研究了漏水量和顶盖压力分布。
分析结果表明水轮机采用适当的止漏环间隙值、装设适当数目的减压排水管和泵板，可以起到减小漏水量、顶盖压力和转轮轴向水推力的效果。
关键词：间隙流道；漏水量；轴向水推力；数值模拟0 引言混流式水轮机的应用最为广泛，这种机组形式具有结构紧凑，运行可靠，能适应很宽的水头范围以及满载时效率高等优点。
但是国内外许多电站在投产运行后出现了不同程度的顶盖压力过高、主轴密封漏水量大和水力振动问题。
工程界对于如何降低顶盖压力和减小主轴密封漏水量，主要是采取设置减压结构的方法。
减压排水管、泵板、引水板和卸荷孔等是比较常用的结构形式。
本文以某高水头电站为研究对象，使用ANSYS FLUENT商业软件，采用RNG模型对密封间隙流道和顶盖上腔部分进行数值模拟，得到止漏环处的漏水量和主轴密封处的压力，计算出减压排水管流速，并分析了密封、减压排水管和泵板的减压效果。
1 数值计算方法本文的计算对象是密封间隙流道和顶盖上腔部分。
采用有限体积法在空间上离散控制方程，在时间离散上使用二阶全隐式格式，压力项应用二阶中心差分格式，其他项采用二阶迎风格式，使用SIMPLEC算法实现压力和速度的分离求解。
边界条件采用压力进口，压力出口条件[2]。
2 物理模型以某高水头电站为研究对象，水轮机的参数是额定水头214.5m，最大水头236.0m，转轮直径1600mm，额定转速600r/min，额定流量16.8 m3/s，额定出力33.2MW。
原始转轮密封间隙为1.35mm，装设2根Ф80mm减压排水管，泵板为径向式，均布8个。
格子玻尔兹曼方法顶盖驱动流
格子玻尔兹曼方法顶盖驱动流格子玻尔兹曼方法（LBM）是一种近年来在流体力学模拟中被广泛应用的数值模拟方法。
该方法可以在复杂的几何和边界条件上精确地模拟各种流体现象。
其中的顶盖驱动流模拟，是一种受到广泛关注的研究领域。
1.什么是顶盖驱动流？顶盖驱动流是指由上端盖板施加外力所形成的流场，它是在一个密闭的正方形边界中进行模拟。
流场通常是由重力和顶部盖板的驱动力组成的。
对于使用LBM计算的顶盖驱动流，由于模型简单，计算效率高，因此已取得了广泛的应用。
2.LBM是如何模拟顶盖驱动流的？LBM的基本原理是通过在大量的离散速度上进行策略性模拟，以模拟物质运动。
与传统的流体动力学方法不同，LBM使用离散的速度和密度来描述流体的运动，因此它可以非常方便的用于复杂的几何和边界条件下的流体模拟。
在顶盖驱动流模拟中，LBM将二维正方形边界分为许多离散化的小单元格，并在每个单元格上施加离散速度和压力值来模拟流体行为。
3.顶盖驱动流模拟的应用领域是什么？顶盖驱动流模拟在许多领域具有广泛的应用，包括地质、生物学、工程学和环境科学等领域。
在地质学中，它被用于模拟岩石的岩石圈运动并对流体流动的地质效应进行分析。
在工程学中，它被用于模拟汽车空气动力学和水力学作用以及结构物的振动和熱传导等现象。
4.顶盖驱动流模拟存在的挑战是什么？尽管顶盖驱动流模拟为模拟流体行为提供了强有力的工具，但仍然存在一些挑战。
例如，该方法需要高度离散化的速度空间和网格结构，这可能会导致计算效率低下和计算成本高昂。
此外，顶盖驱动流模拟还需要对物理设置和计算参数进行大量的调整和测试，才能使模拟结果更加准确和可靠。
总之，格子玻尔兹曼方法顶盖驱动流是一种新兴的数值模拟方法，在流体力学和其他领域的应用越来越广泛。
将来，随着计算机硬件和软件的不断发展，顶盖驱动流模拟将进一步提高计算精度和计算效率，为工程学、生物学和环境科学等领域能提供更准确的解决方案。
格子玻尔兹曼方法顶盖驱动流
格子玻尔兹曼方法顶盖驱动流
格子玻尔兹曼方法是一种用于数值模拟流体动力学的方法，适用于微观尺度上的非平衡流体系统。
其中一个经典的案例是顶盖驱动流。
顶盖驱动流是一种实验室中常见的流体力学现象，其涉及到一个长方形的容器，其中底部被加热，顶部被保持在常温下。
流体在容器内自然对流，并产生了一个横向的流动。
这种流动很适合用格子玻尔兹曼方法来模拟。
在模拟过程中，首先需要将流体分割成小的区域，然后在每个区域中计算流体的速度、密度和分布函数。
通过这些值的计算和更新，可以模拟出流体的运动和演化过程。
针对顶盖驱动流，需要在容器的顶部施加一个固定的速度场，以模拟顶部的运动。
通过格子玻尔兹曼方法的模拟，可以得出流体的速度场、温度场和压力场等信息，从而更好地理解和预测顶盖驱动流的行为。
这对于流体力学的研究和应用都具有重要意义。
- 1 -。
OpenFOAM顶盖驱动流详解!使用手册
OpenFOAM顶盖驱动流详解!使用手册引言这是开源场运算和操作c++库类（openfoam）的使用指南。
他详细描述了OpenFOAM的基本操作。
首先通过第二章一系列教程练习。
然后通过对更多的独立组件的更详细的描述学习openfoam。
Of 首先主要是一个c++库类，主要用于创建可执行文件，比如应用程（application）。
应用程序分成两类：求解器，都是为了解决特定的连续介质力学问题而设计的；公用工程，这些是为了执行包括数据操作等任务而设计的。
Of 包括了数量众多的solver和utilities，牵涉的问题也比较广泛。
将在第三章进行详尽的描述。
Of 的一个强项是用户可以通过必要的预备知识(包括数学，物理和编程技术)创建新的solvers 和utilities。
Of 需要前处理和后处理环境。
前处理、后处理接口就是of本身的实用程序（utilities），以此确保协调的数据传输环境。
图是of总体的结构。
第4章和第五章描述了前处理和运行of 的案例。
既包括用of提供的mesh generator划分网格也包括第三方软件生成的网格数据转换。
第六章介绍后处理。
Chapter 2指导手册在这一章中我们详细描述了安装过程，模拟和后进程处理一些OpenFOAM测试案例，以引导用户运行OpenFOAM的基本程序。
$FOAM_TUTORIALS 目录包含许多案件演示of提供的所有求解器以及许多共用程序的使用，在试图运行教程之前，用户必须首先确保他们已经正确地安装了OpenFOAM。
该教程案件描述 blockMesh预处理工具的使用，paraFoam案例设置和运行OpenFOAM求解器及使用paraFoam进行后处理。
使用OpenFOAM支持的第三方后处理软件的用户可以选择：他们要么可以按照教程使用paraFoam，或当需要后处理时参阅第六章的第三方软件使用说明。
OpenFOAM安装目录下的tutorials目录中所有的指导手册都是可复制的。
方腔顶盖驱动流动
一、问题描述方腔顶盖驱动流动如图1所示的一个简化两维方腔（高，宽都等于L)，内部充满水分。
上表面为移动墙，非维化速度为ｕ／ｕ0 =1。
其他三面为固定墙。
试求方腔内水分流动状态。
u=1, v=0u=0, v=0 u=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 , 称为局部截断误差.截断误差: 实际上，y (x n ) ? y n ， y n 也有误差，它对y n +1的误差也有影响，见下图。
但这里不考虑此误差的影响，仅考虑方法或公式本身带来的误差，因此称为方法误差或截断误差。
局部截断误差的分析：由于假设y n = y (x n ) ，即y n 准确，因此分析局部截断误差时将y (x n +1) 和 y n +1都用点x n 上的信息来表示，工具：Taylor 展开。
显式欧拉公式一阶向前差商近似一阶导数推导如下：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 +++'''=-=+++-+''=+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 +++'≈+↑≈≈=+隐式欧拉公式xn +1点向后差商近似导数 推导如下：几何意义设已知曲线上一点 P n (x n , y n ),过该点作弦线，斜率为(x n +1 , y n +1 ) 点的方向场f (x ,y )方向,若步长h 充分小，可用弦线和垂线x =x n +1的交点近似曲线与垂线的交点。
LB-SGS方法和MRT-LBM方法对高雷诺数顶盖驱动流的数值模拟对比
LB-SGS方法和MRT-LBM方法对高雷诺数顶盖驱动流的数值模拟对比陈春媚;王东;杨志刚【摘要】针对格子Boltzmann方法(Lattice Boltzmann Method,LBM)广泛采用的LBGK模型虽然简单易行,但对高雷诺数流动模拟稳定性不佳的问题,分别采用结合亚格子模型的LBM( LBM with Sub Grid Scale (SGS),LB-SGS)方法和多松驰时间LBM( Multiple Relaxation Time (MRT) LBM,MRTLBM)方法对高雷诺数顶盖驱动流进行数值模拟.取Re =7 500对比2种方法得到的涡位置与标准解之间的误差,结果表明LB-SGS方法更接近标准解；保持雷诺数和顶盖速度不变并减少格点数观察收敛情况,结果表明MRT-LBM方法更稳定.【期刊名称】《计算机辅助工程》【年(卷),期】2012(021)001【总页数】5页(P32-35,49)【关键词】顶盖驱动流;高雷诺数;格子Boltzmann方法;多松驰时间;亚格子模型【作者】陈春媚;王东;杨志刚【作者单位】同济大学上海地面交通工具风洞中心,上海201804;同济大学上海地面交通工具风洞中心,上海201804;同济大学上海地面交通工具风洞中心,上海201804【正文语种】中文【中图分类】O354格子Boltzmann方法(Lattice Boltzmann Method，LBM)源于20世纪70年代提出并发展的格子气自动机方法，可被视为求解连续Boltzmann方程的离散格式之一，又称为格子Boltzmann方程(Lattice Boltzmann Equation，LBE)方法，是流体计算的新方法之一．理论上，一个流体系统可用微观分子动力学、介观动力学模型或宏观连续守恒方程进行描述，同时，也存在一些用这3类方法都不能很好地进行描述的问题．在任一宏观体系中，每个分子的微观运动都遵守力学规律，因此只需计算出大量粒子的个别运动，就可以确定系统的宏观参数，这是分子动力学模拟的基本出发点;另一方面，Boltzmann方程的基本思想为，可以不确定每个分子的运动状态，而是求出每个分子处在某一状态下的概率，通过统计方法得出系统的宏观参数．［1］设速度分布函数为 f，其为空间位置矢量r(x，y，z)，分子速度矢量ξ(ξx，ξy，ξz)以及时间 t的函数．f(r，ξ，t)drdξ 表示 t时刻在 r与 r+dr间的体积元dr=dxdydz中速度在ξ与ξ+dξ间的分子数．分布函数f的控制方程，即为Boltzmann方程式中:a为加速度;Ωf表示碰撞项．QIAN等提出的DdQm(d维空间，m个离散速度)系列模型［1］是 LBM的基本模型．本文采用D2Q9模型(二维九速模型)，其速度配置见图1．eα为图1中9个方向的速度矢量，其中，α=0，1，2，3，4，5，6，7，8，即e0=(0，0)，e1=(1，0)，e2=(0，1)，e3=(－1，0)，e4=(0，－1)，e5=(1，1)，e6=(－1，1)，e7=(－1，－1)，e8=(1，－1)．平衡态分布函数为式中:cs为格子声速;ωα为权因数，顶盖驱动流是计算流体和计算传热学中的经典问题之一，常用作不可压缩流动的校核算例．在顶盖驱动流中，方腔的上边界以某个恒定速度水平右移，其他3个边界保持静止不动．其基本特征为:当流动稳定后，在方腔中央会出现一个一级大涡，而在左下角和右下角会分别出现一个二级涡，当雷诺数Re超过某个临界值后，在方腔左上角还会出现一个涡．这些涡的中心位置是Re的函数，L为方腔的高和宽，U为顶盖的移动速度．顶盖驱动流示意见图2．，若固定网格数和速度不变，随着Re增大黏度减小，流场中的大速度梯度会引起计算的不稳定;另一方面，虽然可以保持黏度和速度，通过增加网格数达到高雷诺数，但过密的网格不但会消耗计算资源，而且容易带来压缩性误差．［2］本文针对高雷诺数，采用结合亚格子模型的LBM(LBM with Sub Grid Scale(SGS)，LBSGS)方法和多松驰时间LBM(Multiple Relaxation Time(MRT)LBM，MRT-LBM)方法对顶盖驱动流进行数值模拟，并比较其准确性和收敛性．在顶盖驱动流中，当Re＞7 500时流动为非定常的，流线随时间变化．由于将涡中心位置作为比较的主要特征，所以本文Re=7 500．演化方程为在LBM的几种模型中，LBGK近似法(又称单松驰因子法)因编程简单、实施容易而备受青睐，但其对边界条件的依赖性强，尤其在高雷诺数时计算稳定性不佳，经常难以收敛．其原因是，LBGK只有一个松驰因子τ=3υ+0．5，其中υ=与LBGK不同的是，式(3)中的松驰因子τ随时间和空间变化，不再是一个常数．根据大涡模拟(Large Eddy Simulation，LES)理论，湍流主要由大涡引起．大涡从流动中获得能量并且各向异性，其之间的相互作用又使能量转化为各向同性的小涡，最后这些小涡将能量耗散成热量．因此在LES中，大涡变量由盒式函数过滤，用非定常N-S方程直接进行数值模拟;采用SGS模型模拟小涡，最常用的 SGS模型为 Smagorinsky模型［3］，由气象学家SMAGORINSKY在1963年提出．LB-SGS方法将经典大涡模拟中的湍流应力模型引入到LBM中，该模型基于涡黏模型理论，采用Smagorinsky模型．由Smagorinsky模型可知式中:Cs＞0，为 Smagorinsky常数，本文取 0．1为流体的变形率张量;τij为亚格子雷诺应力张量;υt为涡黏因数;Δ为滤波宽度，由于本文网格长度统一，过滤宽度即为单位网格长度．LB-SGS模型的运动黏度为松驰因子需要指出的是，在传统方法中变形率张量Sij由差分方程获得;在LBM中，可用分布函数的非平衡项获得，即在 LBGK模型提出后不久，法国学者D’HUMERIERS提出一种广义 LBE 模型(Generalized LBE，GLBE)，该模型曾长期被忽视．2000年，LALLEMAND 和LUO对该模型进行细致的理论分析，表明该模型在物理原理、参数选取和数值稳定性等方面都有很大优势，其与LBGK模型的主要区别在于其碰撞过程使用多个松驰因子，因此被称为 MRT 模型．［4］GLBE可表示为式中:－Λij为碰撞矩阵．D2Q9-MRT模型仍采用规则的正方形格子．变换矩阵为Λ=M－1SM，其中S为一个行数和列数均为9的对角矩阵，对角线元素为(0，se，sε，0，sq，0，sq，sυ，sυ)．剪切黏性和体黏性因数分别为本文所取的松驰因子［5-6］为:se=1．90，sε =LBM的边界处理有很多种，它们起着重要作用，是LBM实施中非常关键的内容之一，会对数值稳定性、计算精度以及计算效率产生很大影响．以下简单介绍本文所用的非平衡外推格式．郭照立等［4］提出一种新的边界处理格式，即非平衡态外推格式，该方法的基本思路为:将边界格点上的分布函数分解为平衡态和非平衡态2部分，根据具体的边界条件构造新的平衡态分布函数，用于近似平衡态部分;对于非平衡态部分，则采用1阶精度的外推方法确定．如图3所示，O为边界格点，B为流体格点，在每次碰撞前需要知道边界格点O上的分布函数fi(O，t)，将其分解为平衡态和非平衡态2部分，即对于平衡态部分(O，t)，如果格点上的宏观物理量已知，则值可以直接得到;如果宏观物理量未知，可由邻近的流体格点B的相应值代替．对于非平衡部分(O，t)，可由流体格点B的非平衡态函数来近似．因此，可得边界格点O的分布函数为非平衡态外推格式主要用于速度边界处理，其整体精度为2阶，具有较好的数值稳定性．正方形方腔高和宽各分布257个格点，顶盖移动速度U=0．1，4条边界均采用非平衡外推法．由图4可知，这2种方法都清楚捕捉到4个涡结构，但从具体位置的坐标上看，LB-SGS更接近一些，数值模拟结果对比见表1．保持雷诺数和顶盖移动速度不变，不断减小空腔的总网格数，观察二者谁先发散．首先采用129×129的网格，在129×129的格点尺寸下，2种方法均捕捉到4个涡结构，见图5．当网格数减少至65×65时，LB-SGS出现发散现象，见图6(a);而MRT-LBM仍能收敛并且清楚捕捉到主涡、左下涡和右下涡，见图6(b)．基于顶盖驱动流，从数值稳定性看，MRT-LBM比LB-SGS更好．(1)针对二维顶盖驱动流问题，分别采用2种LBM模型进行Re=7 500的高雷诺数数值模拟．在2个算例都收敛的情况下，从涡中心位置对比结果可知二者与标准值之间都相差不大，相比之下，LBSGS方法更接近一些;但是在保持雷诺数和顶盖移动速度的情况下，通过减少宽和高的格点数观察二者的收敛情况，MRT-LBM 表现出更好的数值稳定性．(2)由于MRT-LBM有3个可以调节的松驰因子，而 Smagorinsky SGS模型的Cs值也不固定，对于不同的模型，这2个模型的参数都可以更有针对性地进行取值以达到更好的数值模拟结果，更多的分析研究可以围绕这些参数展开．(3)本文的LB-SGS方法采用VC++编程，由于MRT-LBM涉及矩阵运算，因此采用MATLAB编程，后处理软件皆为TECPLOT．MATLAB拥用强大的矩阵运算能力，但执行FOR语句较慢;VC++虽然执行FOR语句快，但涉及矩阵运算需另外编程，二者若能联合仿真，不但可以减少程序的复杂性，而且可以加快程序执行的速度．(4)作为一种新颖的计算流体力学方法，LBM算法简单，程序适合大规模并行计算，可以充分利用多核处理器的资源提升计算性能．随着计算机的发展，该方法将会得到更广泛的应用．王东(1972—)，男，黑龙江哈尔滨人，副教授，博士，研究方向为汽车空气动力学、车辆热管理、气动噪声与计算流体力学，(E-mail)wangdong@tongji．edu．cn【相关文献】［1］何雅玲，王勇，李庆．格子Boltzmann方法的理论及应用［M］．北京:科学出版社，2008:31-55．［2］YANG Fan，LIU Shuhong，WU Yulin，et al．A lattice Boltzmann subgrid model forlid-driven cavity flow［J］．J Hydrodynamics，2005，17(3):289-294．［3］张兆顺，崔桂香，许春晓．湍流大涡数值模拟的理论和应用［M］．北京:清华大学出版社，2008:91-129．［4］郭照立，郑楚光．格子Boltzmann方法的原理及应用［M］．北京:科学出版社，2008:16-24．［5］MUSSA A，ASNARI P，LUO Lishi．Lattice Boltzmann simulations of 2D laminarflows past two tandem cylinders［J］．J Comput Phys，2009，228(4):983-999．［6］YU Huidan，LUO Lishi，GIRIMAJI S S．LES of turbulent square jet flow using an MRT lattice Boltzmann model［J］．Computers＆ Fluids，2006，35(8-9):957-965．。
OpenFOAM顶盖驱动流详解!使用手册(中文翻译版)
OpenFOAM顶盖驱动流详解！使用手册（中文翻译版）引言这是开源场运算和操作 c++库类（openfoam）的使用指南。
他详细描述了OpenFOAM 的基本操作。
首先通过第二章一系列教程练习。
然后通过对更多的独立组件的更详细的描述学习openfoam。
Of 首先主要是一个c++库类，主要用于创建可执行文件，比如应用程序（application）。
应用程序分成两类：求解器，都是为了解决特定的连续介质力学问题而设计的；公用工程，这些是为了执行包括数据操作等任务而设计的。
Of 包括了数量众多的solver和utilities，牵涉的问题也比较广泛。
将在第三章进行详尽的描述。
Of 的一个强项是用户可以通过必要的预备知识(包括数学，物理和编程技术)创建新的solvers 和utilities。
Of 需要前处理和后处理环境。
前处理、后处理接口就是of本身的实用程序（utilities），以此确保协调的数据传输环境。
图1.1是of总体的结构。
第4章和第五章描述了前处理和运行of 的案例。
既包括用of提供的mesh generator划分网格也包括第三方软件生成的网格数据转换。
第六章介绍后处理。
Chapter 2指导手册在这一章中我们详细描述了安装过程，模拟和后进程处理一些OpenFOAM测试案例，以引导用户运行OpenFOAM的基本程序。
$FOAM_TUTORIALS 目录包含许多案件演示of提供的所有求解器以及许多共用程序的使用，在试图运行教程之前，用户必须首先确保他们已经正确地安装了OpenFOAM。
该教程案件描述 blockMesh预处理工具的使用，paraFoam案例设置和运行OpenFOAM 求解器及使用paraFoam进行后处理。
使用OpenFOAM支持的第三方后处理软件的用户可以选择：他们要么可以按照教程使用paraFoam，或当需要后处理时参阅第六章的第三方软件使用说明。
OpenFOAM安装目录下的tutorials目录中所有的指导手册都是可复制的。
- 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。