方腔顶盖驱动流流场数值预测

一、二、问题描述方腔顶盖驱动流动如图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的交点近似曲线与垂线的交点。

格子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。

流场模拟方法流场模拟方法是一种重要的科学技术手段，用于研究和预测流体在各种条件下的运动和相互作用。

它在许多领域中都具有重要应用，如天气预报、风洞试验、环境工程和生物医学研究等。

流体力学是研究流体力学行为的学科，其中流场模拟方法是一个关键的研究领域。

流场模拟方法可以通过数学模型和计算机仿真来预测和分析流体流动的物理特性，从而为各种应用提供有效的解决方案。

流场模拟方法主要包括数值模拟和实验模拟两种。

数值模拟方法是通过建立数学模型和使用计算机算法来模拟流体运动。

这种方法的优点是可以准确预测流场的各种性质，如速度、压力、温度等，并能够在很短的时间内得到结果。

然而，数值模拟方法需要依赖复杂的数学模型和计算机算法，因此对计算资源要求高，而且模拟结果可能受到模型的假设和参数选择的影响。

实验模拟方法是通过设计和进行实验来模拟流体运动。

这种方法的优点是可以直接观测和测量流体的运动和相互作用，对结果的可信度高。

同时，实验模拟方法也能够提供丰富的数据来验证和改进数值模拟方法。

然而，实验模拟方法需要大量的设备和实验操作，并且受到实验条件和测量误差的限制。

在流场模拟方法中，数值模拟方法常用的技术包括有限元法、有限差分法和边界元法等。

这些技术通过对流体运动的偏微分方程进行离散化和求解，从而获得流场的数值解。

有限元法是一种广泛应用的数值模拟方法，它把流场划分为多个小单元，然后通过求解各单元上的方程来获得整个流场的数值解。

有限差分法是另一种常用的数值模拟方法，它将流场划分为网格点，在每个网格点上计算流体的变化量，然后通过迭代求解来获得整个流场的数值解。

边界元法是一种基于边界条件的数值模拟方法，它将流场划分为多个边界元，然后通过求解边界元上的方程来获得整个流场的数值解。

这些数值模拟方法都有各自的优点和适用范围，在具体应用中需要根据问题的复杂程度和计算资源的限制来选择合适的方法。

实验模拟方法中常用的技术包括风洞试验、流体力学实验和粒子图像测速法（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 顶盖驱动方腔流动数值模拟研究作者：杨约翰
来源：《科学与财富》2020年第20期
摘要：稀薄气体条件下，气体分子的宏观性质会受到气体分子的间断粒子效应的影响，连续介质模型不能完全反映物理实际。

建立在连续介质模型基础上的流体力學控制方程 NS 方程就不能适用于此。

本文选用 EGKS 模型的 Blotzmann 方程作为其控制方程进行了数值模拟，对不同 Kn 数下顶盖驱动方腔流动的流动特性进行了研究。

关键词：稀薄气体;NS 方程;Blotzmann 方程;EGKS 模型;顶盖驱动方腔流动。

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）是一种常用于表征复杂流体动力学行为的流体力学方法。

顶盖驱动流谱方法通过在流体上方施加一个顶盖运动，并通过监测流体中的扩散和剪切变形来评估流体的流动性质。

顶盖驱动流谱方法的基本原理是利用顶盖的运动来引起流体的流动。

通过在流体上方施加一个顶盖的运动，可以在流体中引入额外的扩散和剪切变形。

顶盖的运动可以通过不同的方式实现，如旋转、振荡或往复移动。

同时，在顶盖运动的影响下，流体中的粒子或标记物会发生位移和旋转，这些运动可以通过光学、电学或其他检测方法进行监测和分析。

顶盖驱动流谱方法可以用于研究不同流体性质的变化，如黏度、流变特性和颗粒浓度等。

通过测量扩散和剪切变形的程度，可以得到关于流体动力学行为的定量信息。

这种方法在生命科学领域的应用也很广泛，可以用于研究细胞运动、颗粒流和微流体中的流动等。

与传统的流体力学方法相比，顶盖驱动流谱方法具有一些优势。

首先，它可以对复杂、非牛顿流体的性质进行研究，而传统方法往往只适用于牛顿流体。

其次，顶盖驱动流谱方法可以用于研究微观尺度下的流体行为，可以探索纳米和微米尺度上的流动现象。

此外，这种方法还可以与其他技术相结合，如荧光标记、微流控等，进一步拓展其应用范围。

总之，顶盖驱动流谱方法是一种有效的流体力学方法，可以用
于研究复杂流体的动力学行为。

其原理简单，应用范围广泛，有望在生物医学和纳米技术等领域发挥重要作用。