Fortran编写的二维不可压定常空腔流动数值程序数值发散
问题遇到的现象和发生背景
我在用Fortran编写模拟二维不可压定常空腔流动的时候一直数值发散，各个网格点的修正压力一直增大，最后趋于无穷了。
流动条件：二维不可压定常流动（Re=100 驱动速度为1，指定了上层壁面为驱动的壁面）
控制方程：无量纲的动量与连续方程
时间推进方法：SMAC
网格：普通的交错网格，在边界处为边界条件补充ghost grid（速度的ghost grid为了空间离散，压力ghost
grid 为了泊松方程）
边界条件：Neumann 边界条件，壁面上的速度二阶导数以及压力一阶导数为0，对壁面采用中心差得到
ghost grid的速度与修正压力信息。
没有质量流进流出再结合壁面无滑移条件与上层的驱动速度，就可以得到各个壁面的速度值的条件了。
泊松方程解法：Gauss迭代并且采用超松弛，现在设定的松弛因子为1.7.迭代步数也进行了限制现在设定为
50000步，但是不知道这个是不是适合。
泊松迭代收敛判断条件是
主循环的收敛条件是根据是否满足连续性方程来判断的。
问题相关代码
接下来是的代码
Program main
implicit none !
Real::control,step,control_step,sor
integer::i,j !坐标
Real::a,b,c !泊松方程系数
integer::nx=32,ny=32 !x与y方向格子数目
Real::re !雷诺数
real::n !计算用时
Real::bctop !上边界驱动速度无量纲
Real::dt !时间步长
Real::dx,dy !x与y方向网格宽度，因为是均匀的，所以不用数组
Real::sum,sum1,sum2
Real,ALLOCATABLE,Dimension(:,:)::u,v,p !u,v,p 可分配数组
Real,ALLOCATABLE,Dimension(:,:)::up,vp !u,v预测速度
Real,ALLOCATABLE,Dimension(:,:)::fa !修正压力
Real,ALLOCATABLE,Dimension(:,:)::dpx,dpy,ddx,ddy,dcx,dcy! 遗迹压力扩散对流项
Real,ALLOCATABLE,Dimension(:,:)::d! 泊松方程系数
ALLOCATE(u(0:nx+1,0:ny+1)) !原本为x:0.5,nx+0.5 y:0,ny+1
ALLOCATE(up(0:nx+1,0:ny+1)) !原本为x:0.5,nx+0.5 y:0,ny+1
ALLOCATE(v(0:nx+1,0:ny+1)) !原本为x:0,nx+1 y:0.5,ny+0.5
ALLOCATE(vp(0:nx+1,0:ny+1)) !原本为x:0,nx+1 y:0.5,ny+0.5
ALLOCATE(p(0:nx+1,0:ny+1)) ! calculation area (1:nx,1:ny)
ALLOCATE(fa(0:nx+1,0:ny+1))
ALLOCATE(dpx(0:nx+1,0:ny+1)) !x:1+0.5到nx-0.5 y:1到ny
ALLOCATE(ddx(0:nx+1,0:ny+1)) !x:1+0.5到nx-0.5 y:1到ny
ALLOCATE(dcx(0:nx+1,0:ny+1)) !x:1+0.5到nx-0.5 y:1到ny
ALLOCATE(dpy(0:nx+1,0:ny+1)) !x:1到nx y:1+0.5到ny-0.5
ALLOCATE(ddy(0:nx+1,0:ny+1)) !x:1到nx y:1+0.5到ny-0.5
ALLOCATE(dcy(0:nx+1,0:ny+1)) !x:1到nx y:1+0.5到ny-0.5
ALLOCATE(d(0:nx+1,0:ny+1)) !x:1到nx y:1+0.5到ny-0.5
! initialize
u=0
v=0
up=0
vp=0
p=0
fa=0
dpx=0
dpy=0
ddx=0
ddy=0
dcx=0
dcy=0
sum=0
sum1=0
sum2=0
sum3=0
fa=0
step=1
control=0
d=0
i=1
j=1
n=1
!setting
dx=1./REAL(nx)
dy=1./REAL(ny)
a=1./(dx**2) !泊松方程系数
c=1./(dy**2) !泊松方程系数
b=2.*a+2.*c !泊松方程系数
sor=1.7 !松弛因子
dt=0.01 !单个时间步长
bctop=1 !上bc
re=100 !雷诺数
control_step=5000000 !迭代步数控制
110 FORMAT(34(F6.2,X))
100 FORMAT(34(F7.2,2X))
Do !主循环
!BC condition
ghostleftright:Do j=1,ny-1 !Neumann BC 左右边界上速度条件 ghost grid
v(0,j)=-v(1,j)
u(0,j)=0
v(nx+1,j)=-v(nx,j)
u(nx+1,j)=0
end do ghostleftright
ghostbottom:Do i=1,nx-1 !Neumann BC 下ghost grid
u(i,0)=-u(i,1)
v(i,0)=0
end do ghostbottom
ghostabove:Do i=0,nx !Neumann BC 上ghost grid
u(i,ny+1)=2.*bctop-u(i,ny)
v(i,ny+1)=0
end do ghostabove
!计算x压力扩散对流
calx:Do i=1,nx-1
do j=1,ny
dpx(i,j)=(p(i+1,j)-p(i,j))/dx
ddx(i,j)=(u(i-1,j)-2*u(i,j)+u(i+1,j))/(dx**2)+(u(i,j-1)-2*u(i,j)+u(i,j+1))/(dy**2)
dcx(i,j)=(((u(i+1,j)+u(i,j))/2.)**2-((u(i,j)+u(i-1,j))/2.)**2)/dx+&
(((u(i,j)+u(i,j+1))/2.)*((v(i,j)+v(i+1,j))/2.)-&
((u(i,j)+u(i,j-1))/2.)*((v(i,j-1)+v(i+1,j-1))/2.))/dy
end do
end do calx
!算u预测速度
Do i=1,nx-1
do j=1,ny
up(i,j)=u(i,j)+dt*(-dcx(i,j)-dpx(i,j)+(ddx(i,j)/re))
End do
End do
!计算y压力扩散对流
caly:Do i=1,nx
do j=1,ny-1
dpy(i,j)=(p(i,j+1)-p(i,j))/dy
ddy(i,j)=(v(i-1,j)-2.*v(i,j)+v(i+1,j))/(dx**2)+(v(i,j-1)-2.*v(i,j)+v(i,j+1))/(dy**2)
dcy(i,j)=(((v(i,j+1)+v(i,j))/2.)**2-((v(i,j)+v(i,j-1))/2.)**2)/dy+&
(((u(i,j)+u(i,j+1))/2.)*((v(i,j)+v(i+1,j))/2.)-&
((u(i-1,j)+u(i-1,j+1))/2.)*((v(i,j)+v(i-1,j))/2.))/dx
end do
end do caly
!算v预测速度
Do i=1,nx
do j=1,ny-1
vp(i,j)=v(i,j)+dt*(-dcy(i,j)-dpy(i,j)+(ddy(i,j)/re))
End do
End do
!计算预测速度差分
do i=1,nx !x方向推进
do j=1,ny !y方向推进
d(i,j)=((up(i,j)-up(i-1,j))/dx+(vp(i,j)-vp(i,j-1)/dy))/dt
End do
End do
!泊松
main_poisson:Do !主迭代循环
Do j=1,ny !Neumann BC 压力梯度为0边界上用中心差分得 ghost grid left right
fa(0,j)=fa(1,j)
fa(nx+1,j)=fa(nx,j)
end do
Do i=1,nx !Neumann BC 压力梯度为0边界上用中心差分得 ghost grid bottom top
fa(i,0)=fa(i,1)
fa(i,ny+1)=fa(i,ny)
end do
outer:do i=1,nx !x方向推进
inner:do j=1,ny !y方向推进
fa(i,j)=fa(i,j)&
+sor*(a*fa(i-1,j)+a*fa(i+1,j)-b*fa(i,j)+c*fa(i,j-1)+c*fa(i,j+1)-d(i,j))/b
End do inner
End do outer
! 计算误差
Do i=1,nx
Do j=1,ny
sum=sum+((a*fa(i-1,j)+a*fa(i+1,j)-b*fa(i,j)+c*fa(i,j-1)+c*fa(i,j+1))&
-d(i,j))**2
sum1=sum1+d(i,j)**2
end do
end do
sum=sqrt(sum/sum1)
if(sum<=1E-5)then !符合收敛要求
control=1
else
step=step+1 !增加步数计数
end if
if (step==control_step)then !到达上限步数，输入err错误
control=1
else
step=step+1
end if
if(control==1)exit !控制主循环
sum=0 !初始化求和
sum1=0
End do main_poisson
write(*,*)step,sum,sum1
Do i=1,nx-1 !u时间推进
do j=1,ny
u(i,j)=up(i,j)-(fa(i+1,j)-fa(i,j))/dx*dt
End do
end do
Do i=1,nx !v时间推进
do j=1,ny-1
v(i,j)=vp(i,j)-(fa(i,j+1)-fa(i,j))/dy*dt
End do
end do
Do i=1,nx !p时间推进
do j=1,ny
p(i,j)=p(i,j)+fa(i,j)
End do
end do
! 以是否满足连续性方程为收敛判据
! u差分
Do i=1,nx
do j=1,ny
sum2=sum2+((u(i,j)-u(i-1,j))/dx+((v(i,j)-v(i,j-1))/dy))**2
end do
end do
write(*,*)"/................../"
write(*,*)"up"
write(*,110)((up(i,j),i=0,nx+1),j=ny+1,0,-1)
write(*,*)"/................../"
write(*,*)"/................../"
write(*,*)"vp"
write(*,110)((vp(i,j),i=0,nx+1),j=ny+1,0,-1)
write(*,*)"/................../"
write(*,*)"/................../"
write(*,*)"u"
write(*,110)((u(i,j),i=0,nx+1),j=ny+1,0,-1)
write(*,*)"/................../"
write(*,*)"/................../"
write(*,*)"v"
write(*,110)((v(i,j),i=0,nx+1),j=ny+1,0,-1)
write(*,*)"/................../"
write(*,*)"p"
write(*,110)((p(i,j),i=0,nx+1),j=ny+1,0,-1)
write(*,*)"/................../"
write(*,*)"fa"
write(*,110)((fa(i,j),i=0,nx+1),j=ny+1,0,-1)
write(*,*)"/................../"
sum2=sqrt(sum2/(nx*ny))
write(*,*)sum4,n
If(abs(sum2)<=1E-5.OR.n==200.OR.err==1)exit
n=n+1 !计算用了几个时间步长
sum=0
sum1=0
sum2=0
up=0
vp=0
d=0
fa=0
temp=0
step=1
control=0
d=0
End do
END PROGRAM
运行结果
修正压力会一直变大，第六步左右就无穷大了。
尝试过的方法
起初我认为是我的语句有问题，但是我检查了好几遍并没有发现错误，之后又对理论进行了检查，还是没有发现问题，之后我查询了下Gauss迭代是否适合求解这类问题，根据我查询到的内容，Gauss迭代是可以求解的。我也更改了各种参数，但是问题还是未能解决。希望各位大大能帮我看下我是在哪个部分有问题。