利用算法和压力迭代解法求解二维不可压缩黏性流动问题Fortran语言

Fortran77语言源程序!MAC-Chorin2D.for/* *利用MAC算法和Chorin压力迭代解法求解二维不可压缩黏性Poisuille流动问题(Fortran77语言版本) * */programMAC_Chorinimplicitdoubleprecision(a-h,o-z)parameter(M1=101,M2=21)common/G_def/Jx,Jy,dx,dy,dt,Re,EPSdimensionu(0:M1+1,0:M2+1),v(0:M1+1,0:M2+1),p(0:M1+1,0:M2+1)dimensiontu(0:M1+1,0:M2+1),tv(0:M1+1,0:M2+1)!x方向网格点数Jx=M1!y方向网格点数Jy=M2!Reynolds数Re=100!x方向网格间距dx=1.0/(Jx-1)!y方向网格间距dy=0.2/(Jy-1)!时间步长dt=0.0002!收敛限EPS=1d-6callinit(u,v,p)!lambda常数与时间步长、网格间距和Reynolds数都有关。关系到收敛性，应仔细选择!dlambda=2d-3;n=0!计算压力1callsolvep(p,u,v,dlambda)!计算速度callsolveuv(u,v,p,tu,tv)!判断是否满足连续性条件isconv=iconv(u,v,dm)n=n+1if(mod(n,100).eq.0)write(*,*)n,'StepMaxVelDivergence=',dmif(isconv.eq.0)goto1calloutput(u,v,p)end! !速度边界条件!入口：无；!出口：u、v，已经给定边界，u、v上下边界为无滑移固壁，u、v左右边界为!出口边界。!subroutinebounduv(u,v)implicitdoubleprecision(a-h,o-z)common/G_def/Jx,Jy,dx,dy,dt,Re,EPSdimensionu(0:Jx+1,0:Jy+1),v(0:Jx+1,0:Jy+1)do10i=0,Jx+1u(i,1)=0;u(i,Jy)=0;v(i,0)=-v(i,1);v(i,Jy)=-v(i,Jy-1);continuedo11j=1,Jyu(0,j)=u(1,j);u(Jx,j)=u(Jx-1,j);v(1,j)=v(1,j);v(Jx,j)=v(Jx-1,j);continueend! !压力边界条件!入口：无；!出口：p,为已经给定边界，只需要给定左右压力，左侧1.0,右侧0.0。! subroutineboundp(p)implicitdoubleprecision(a-h,o-z)common/G_def/Jx,Jy,dx,dy,dt,Re,EPSdimensionp(0:Jx+1,0:Jy+1)do20j=0,Jy+1p(1,j)=1.0;p(Jx,j)=0.0;20continueend! !初始化!入口：无；!出口：u、v、p，已经给定的初始值。! subroutineinit(u,v,p)implicitdoubleprecision(a-h,o-z)common/G_def/Jx,Jy,dx,dy,dt,Re,EPSdimensionu(0:Jx+1,0:Jy+1),v(0:Jx+1,0:Jy+1),p(0:Jx+1,0:Jy+1)do30i=0,Jx+1do30j=0,Jy+1u(i,j)=0;v(i,j)=0;p(i,j)=0;30continuecallbounduv(u,v)callboundp(p);end! !根据n时刻u、v、p.采用中心差分格式计算n+1时刻u、v!入口：u、v、p,n时刻u、v、p,tu、tv临时变量；!出口：u、v,n+1时刻u、v。! subroutinesolveuv(u,v,p,tu,tv)implicitdoubleprecision(a-h,o-z)common/G_def/Jx,Jy,dx,dy,dt,Re,EPSdimensionu(0:Jx+1,0:Jy+1),v(0:Jx+1,0:Jy+1),p(0:Jx+1,0:Jy+1)dimensiontu(0:Jx+1,0:Jy+1),tv(0:Jx+1,0:Jy+1)do40i=1,Jx-1do40j=2,Jy-1vav=(v(i,j)+v(i+1,j)+v(i,j-1)+v(i+1,j-1))/4.0;adv=-u(i,j)*(u(i+1,j)-u(i-1,j))/dx/2.0$ -vav*(u(i,j+1)-u(i,j-1))/dy/2.0;vis=((u(i+1,j)-2*u(i,j)+u(i-1,j))/dx/dx$ +(u(i,j+1)-2*u(i,j)+u(i,j-1))/dy/dy)/Re;prs=-(p(i+1,j)-p(i,j))/dx;tu(i,j)=(adv+vis+prs)*dt+u(i,j);continuedo41i=2,Jx-1do41j=1,Jy-1uav=(u(i-1,j+1)+u(i,j+1)+u(i-1,j)+u(i,j))/4.0;adv=-uav*(v(i+1,j)-v(i-1,j))/dx/2.0$ -v(i,j)*(v(i,j+1)-v(i,j-1))/dy/2.0;vis=((v(i+1,j)-2*v(i,j)+v(i-1,j))/dx/dx$ +(v(i,j+1)-2*v(i,j)+v(i,j-1))/dy/dy)/Re;prs=-(p(i,j+1)-p(i,j))/dy;tv(i,j)=(adv+vis+prs)*dt+v(i,j);continuedo42i=1,Jxdo42j=1,Jyu(i，j)=tu(i,j);v(i,j)=tv(i,j);continuecallbounduv(u,v)end! !根据n时刻u、v、p计算n+1时刻p!入口：u、v、p,n时刻u、v、p；dlambda，常数；!出口:p,n+1时刻p。! subroutinesolvep(p,u,v,dlambda)implicitdoubleprecision(a-h,o-z)common/G_def/Jx,Jy,dx,dy,dt,Re,EPSdimensionu(0:Jx+1,0:Jy+1),v(0:Jx+1,0:Jy+1),p(0:Jx+1,0:Jy+1)do50i=1,Jxdo50j=1,Jyp(i,j)=p(i,j)$ -dlambda*((u(i,j)-u(i-1,j))/dx+(v(i,j)-v(i,j-1))/dy);50continuecallboundp(p)end! !判断收敛!入口：u、v,当前时刻u、v、p；!出口：dm,当前网格点上速度散度最大值，!返回：1达到收敛，0尚未收敛。! integerfunctioniconv(u,v,dm)implicitdoubleprecision(a-h,o-z)common/G_def/Jx,Jy,dx,dy,dt,Re,EPSdimensionu(0:Jx+1,0:Jy+1),v(0:Jx+1,0:Jy+1),p(0:Jx+1,0:Jy+1)dm=0;do60i=1,Jxdo60j=1,Jyer=dabs((u(i,j)-u(i-1,j))/dx+(v(i,j)-v(i,j-1))/dy)if(dm.lt.er)dm=er60continueif(dm.lt.EPS)theniconv=1elseiconv=0endifend! !输出全场计算结果和x=0.5处的速度u剖面!入口：u、v、p，需要输出的u、v、p；!出口：无。! subroutineoutput(u,v,p)implicitdoubleprecision(a-h,o-z)common/G_def/Jx,Jy,dx,dy,dt,Re,EPSdimensionu(0:Jx+1,0:Jy+1),v(0:Jx+1,0:Jy+1),p(0:Jx+1,0:Jy+1)open(1,file='result.plt',status='unknown')write(1,*)'Title="computationalResults"'write(1,*)'VARIABLES="x","y","u","v","p"$ZONET="Zone1"'write(1,*)'I=',Jy,'J=',Jx,'K=',1,'ZONETYPE=Ordered'write(1,*)'DATAPACKING=POINT'do70i=1,Jxdo70j=1,Jywrite(1,71)(i-1)*dx,(j-1)*dy,(u(i-1,j)+u(i,j))/2.0,$ (v(i,j-1)+v(i,j))/2.0,p(i,j)70continueclose(