以一階迎風(fēng)型離散格式和壓力迭代求解二維流動問題Fortran語言

1、Fortran77語言源程序program fvm_upwind_MAC_couette!以一階迎風(fēng)型離散格式和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.!雷諾數(shù)dt=0.0005!時間步長lambda=0.002 !Chor

2、in壓力迭代常數(shù)，直接影響收斂性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! 一階迎風(fēng)型離散格式!輸入：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數(shù)據(jù)格式畫圖!輸入：u、v、p、mx、my、dx、dy，最后結(jié)果；!輸出：無。!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