以一阶迎风型离散格式和压力迭代求解二维流动问题Fortran语言

1、Fortran77语言源程序program fvm_upwind_MAC_couette!以一阶迎风型离散格式和Chorin压力迭代求解二维Couette流动问题(Fortran77语言版本)!parameter(mx=101,my=21)implicit double precision(a-h,o-z)dimension u(mx,my+1),v(mx+1,my),p(mx+1,my+1),un(mx,my+1),vn(mx+1,my)double precision lambdainteger stepre=100.!雷诺数dt=0.0005!时间步长lambda=0.002 !Chor

2、in压力迭代常数，直接影响收敛性dx=1.0d0/(mx-1)dy=0.2d0/(my-1)step=0err=100.call initial(u,v,p)!初始化do while(err1e-6.and.steperr) err=tempu(i,j)=un(i,j)enddoenddodo i=1,mx+1do j=1,mytemp=abs(v(i,j)-vn(i,j)/dtif(temperr) err=tempv(i,j)=vn(i,j)enddoenddowrite(*,*) err=,errenddocall output(u,v,p,mx,my,dx,dy) !输出End!初始化

3、!输入：无；!输出：u、v、p,初始速度、压强。!subroutine initial(u,v,p)parameter(mx=101,my=21)double precision u(mx,my+1),v(mx+1,my),p(mx+1,my+1)do i=1,mx+1do j=1,my+1p(i,j)=1.0enddoenddodo i=1,mxdo j=1,my+1u(i,j)=0.0if(j=my+1) u(i,j)=2.0enddoenddodo i=1,mx+1do j=1,myv(i,j)=0.0enddoenddoendsubroutine! Chorin 压力迭代!输入：u、v

4、、p、mx、my、dx、dy、lambda , n时刻速、压强、其他各量;!输出：p，n+1时刻压强。!subroutine calp(u,v,p,mx,my,dx,dy,lambda)implicit double precision(a-h,o-z)dimension u(mx,my+1),v(mx+1,my),p(mx+1,my+1),d(mx+1,my+1)double precision lambdado i=2,mxdo j=2,myd(i,j)=(u(i,j)-u(i-1,j)/dx+(v(i,j)-v(i,j-1)/dyenddoenddodo i=2,mxdo j=2,my

5、p(i,j)=p(i,j)-lambda*d(i,j)enddoenddodo i=2,mxp(i,1)=p(i,2)p(i,my+1)=p(i,my)enddodo j=1,my+1p(1,j)=p(2,j)p(mx+1,j)=p(mx,j)enddoendsubroutine! 一阶迎风型离散格式!输入：u、v、pn、mx、my、dx、dy、dt、re , n时刻速度、n+1时刻压强、其他各量；!输出：un,vn，n+1时刻速度。!subroutine upwind(u,v,pn,un,vn,mx,my,dx,dy,dt,re)implicit double precision(a-h,o

6、-z)dimension u(mx,my+1),v(mx+1,my),pn(mx+1,my+1),un(mx,my+1),vn(mx+1,my)double precision miumiu=1.0/redo i=2,mx-1do j=2,myaw=miu+max(0.5*(u(i-1,j)+u(i,j)*dy,0.0)ae=miu+max(0.0,-0.5*(u(i,j)+u(i+1,j)*dy) as=miu+max(0.5*(v(i,j-1)+v(i+1,j-1)*dx,0.0) an=miu+max(0.0,-0.5*(v(i,j)+v(i+1,j)*dx) df=0.5*(u(i+1

7、,j)-u(i-1,j)*dy+0.5*(v(i,j)+v(i+1,j)-v(i,j-1)-v(i+1,j-1)*dx ap=aw+ae+as+an+df un(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) & -dt*(pn(i+1,j)-pn(i,j)/dxenddoenddo!do i=2,mx-1un(i,1)=-un(i,2)un(i,my+1)=2.0-un(i,my)enddodo j=1,my+1un(1,j)=un(2,j)un(mx,j)=un(mx-1,j)

8、enddo!u边界条件!do i=2,mxdo j=2,my-1aw=miu+max(0.5*(u(i-1,j)+u(i-1,j+1)*dy,0.0)ae=miu+max(0.0,-0.5*(u(i,j)+u(i,j+1)*dy) as=miu+max(0.5*(v(i,j-1)+v(i,j)*dx,0.0) an=miu+max(0.0,-0.5*(v(i,j)+v(i,j+1)*dx) df=0.5*(u(i,j)+u(i,j+1)-u(i-1,j)-u(i-1,j+1)*dy+0.5*(v(i,j+1)-v(i,j-1)*dx ap=aw+ae+as+an+dfvn(i,j)=v(i,j

9、)+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) & -dt*(pn(i,j+1)-pn(i,j)/dyenddoenddo!do i=2,mxvn(i,1)=0.0vn(i,my)=0.0enddodo j=1,myvn(1,j)=-vn(2,j)vn(mx+1,j)=-vn(mx,j)enddo!v边界条件!Endsubroutine!格式输出，采用Tecplot数据格式画图!输入：u、v、p、mx、my、dx、dy，最后结果；!输出：无。!subroutine output(u,v,p,mx,my

10、,dx,dy)implicit double precision(a-h,o-z)dimension u(mx,my+1),v(mx+1,my),p(mx+1,my+1)dimension uc(mx,my),vc(mx,my),pc(mx,my),x(mx),y(my)do i=1,mxx(i)=(i-1)*dxenddodo j=1,myy(j)=(j-1)*dyenddodo i=1,mxdo j=1,myuc(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+1,j)+p(i,j+1)+p(i+1,j+1)enddoenddoopen(1,file=result.plt)write(1,*) TITLE = 2D Resultswrite(1,*) VARI