dvm_separate_point.f90

圆柱绕流问题的涡方法数值模拟（fortran源码）

subroutine seprate_point()
!计算KNP（剥离点数目）以及剥离点几何位置
use varible
integer ::i,j,k,m,l,self
real ::spv(knp,2)
real dis

!-------------------------剥离点数目及编号计算--------------------------
entry nsp()
nt(1)=1
nt(2)=5
nt(3)=10
nt(4)=31
nt(5)=36
do i=1,knp
if(nt(i)/=0)then
num=nt(i)
bnp(num)=1
end if
end do

return

entry spini()
!---------------------剥离点未剥离时的几何参数----------------------
  do i=1,knp
   k=nt(i)
   tao(1,i)=btao(k)
   Gx(1,i)=bv(k,1)*(1+pi/nall)
   Gy(1,i)=bv(k,2)*(1+pi/nall)
  end do
!---------------------剥离点一个dt时间步长内的初始计算
  do m=1,mpv
     do i=1,knp
	 self=nt(i)
	 spv(i,1)=0.0
	 spv(i,2)=0.0
	    do j=1,nall
	       if(self==j)then 
	         continue
	       else
		     dis=(Gx(1,i)-bv(j,1))**2+(Gy(1,i)-bv(j,2))**2
			 spv(i,1)=spv(i,1)+btao(j)*(-(Gy(1,i)-bv(j,2)))/2.0/pi/dis
		     spv(i,2)=spv(i,2)+btao(j)*(Gx(1,i)-bv(j,1))/2.0/pi/dis
		   end if
        end do
           spv(i,1)=U*cos(alpha)-spv(i,1)
		   spv(i,2)=U*sin(alpha)-spv(i,2)	    
	 end do

	 do i=1,knp
	 Gx(1,i)=Gx(1,i)+spv(i,1)*ddt
	 Gy(1,i)=Gy(1,i)+spv(i,2)*ddt
	 end do
! open(unit=10,file='GxGy.txt')
	 !write(*,*)(gx(1,i),i=1,knp)
	 !write(*,*)(gy(1,i),i=1,knp)
  end do
!  write(*,*)(gx(1,i),i=1,knp)
return	   
end subroutine