backup.f90

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

 do k=1,kk
      do L=1,knp
      dx(k,L)=Gx(k,L)
      dy(k,L)=Gy(k,L)
      end do
   end do
   do i=1,nall
   fa(i)=btao(i)
   end do
!--------------------壁面镜相-------------------

!-------------------新生边界涡强度计算-------------------
   do i=1,nall
   tempv=0
      do L=1,knp
	     do k=1,kk
		 if (ndp(k,L)==0)then
		  continue 
		 else
		  dis=sqrt((bc(i,1)-dx(k,L))**2+(bc(i,2)-dy(k,L))**2)
		  if(dis<sigma(k))then 
		   dis=sigma(k)
		  else
		   continue
		  end if
		  tempv=tempv+tao(k,L)*(psym(i,1)*(bc(i,1)-dx(k,L))+psym(i,2)*(bc(i,2)-dy(k,L)))/2.0/pi/dis/dis
		 end if
		 end do
      end do
   bphy(i)=hsym(i,1)*u*cos(alpha)+hsym(i,2)*u*sin(alpha)-tempv
   end do
   bphy(nmax)=0.0
   do L=1,knp
      do k=1,kk
	  bphy(nmax)=bphy(nmax)+tao(k,L)
	  end do
   end do
   bphy(nmax)=-bphy(nmax)
!--------------去掉一个多余条件---------------
   do i=1,nall
      do j=1,nall 
      bmatrx(i,j)=0
	     do k=1,nmax
		 bmatrx(i,j)=bmatrx(i,j)+matrx1(k,i)*matrx1(k,j)
		 end do
	  end do
   end do
   do i=1,nall
   btao(i)=0
      do k=1,nmax
      btao(i)=btao(i)+matrx1(k,i)*bphy(k)
	  end do
   end do
!------------计算新一轮壁面的壁面涡强度----------------------
   call gs(bmatrx,btao,nall)
!   write(*,*)(btao(i),i=1,nall)
!------------计算KK时刻所有涡之间的诱导速度及点涡新坐标------------------

   do m=1,mpv
      do L=1,knp
	     do k=1,kk
		    VAx(k,L)=0.0
		    VAy(k,L)=0.0
            if (ndp(k,L)==0)then
		      continue
		    else
		      do i=1,nall
		      dis=sqrt((dx(k,L)-bv(i,1))**2+(dy(k,L)-bv(i,2))**2)
		      if (dis<sigma0(k,M))then
		        dis=sigma0(k,M)
	   	      else 
		        continue
		      end if
		      VAx(k,L)=VAx(k,L)+btao(i)*(-(dy(k,L)-bv(i,2)))/2/pi/dis/dis
		      VAy(k,L)=Vay(k,L)+btao(i)*(dx(k,L)-bv(i,1))/2/pi/dis/dis
		      end do
		    end if

		    do k1=1,kk
		    if (k1==k)then
	          continue
			else 
			  if (ndp(k1,L)==0)then
			    continue
			  else
			    dis=sqrt((dx(k,L)-dx(k1,L))**2+(dy(k,L)-dy(k1,L))**2)
			    if (dis<sigma0(k1,M))then
			      dis=sigma0(k1,M)
			    else
			      continue 
			    end if
			    VAx(k,L)=VAx(k,L)+tao(k1,L)*(-(dy(k,L)-dy(k1,L)))/2/pi/dis/dis
				VAy(k,L)=VAx(k,L)+tao(k1,L)*(dx(k,L)-dx(k1,L))/2/pi/dis/dis
			  end if
			end if
			end do

			do L1=1,Knp
			if (L1==L)then
			  continue
			else
			  do k2=1,kk
			  if (ndp(k2,L1)==0)then
			    continue
			  else
			    dis=sqrt((dx(k,L)-dx(k2,L1))**2+(dy(k,L)-dy(k2,L1))**2)
			    if (dis<sigma0(k2,M))then
				  dis=sigma0(k2,M)
				else 
				  continue
				end if
				VAx(k,L)=VAx(k,L)+tao(k2,L1)*(-(dy(k,L)-dy(k2,L1)))/2/pi/dis/dis
				VAy(k,L)=VAy(k,L)+tao(k2,L1)*(dx(k,L)-dx(k2,L1))/2/pi/dis/dis
		      end if
			  end do
			end if
			end do

			VAx(k,L)=U*cos(alpha)-VAx(k,L)
			VAy(k,L)=U*sin(alpha)-VAx(k,L)
		 end do
	  end do
	  do L=1,knp
	     do k=1,kk
		 if (ndp(k,L)==0)then 
	       continue
		 else
		   dx(k,L)=dx(k,L)+VAx(k,L)*ddt
		   dy(k,L)=dx(k,L)+VAx(k,L)*ddt
		 end if 
		 end do
	  end do
	  call combination(kk,dx,dy)
   end do


!--------------------壁面镜相-------------------
do L=1,knp
   do k=1,kk
   if(ndp(k,L)==0)then
   continue
   else
     if (dx(k,L)<xl2.and.abs(dy(k,L))<yl2)then
	 dx(k,L)=(xl2-dx(k,L))+xl2
	 else
	 continue
	 endif
   end if
   end do
end do
!   do k=1,kk
!      do L=1,knp
!	     if((xl2-0.02)<=dx(k,L).and.dx(k,L)<=xl2.and.(-yl2+0.02)<=dy(k,L).and.dy(k,L)<=(yl2-0.02))then
!		   dx(k,L)=xl2-dx(k,L)+xl2
!		 elseif(-xl2<=dx(k,L).and.dx(k,L)<=(-xl2+0.02).and.(-yl2+0.02)<=dy(k,L).and.dy(k,L)<=(yl2-0.02))then
!		   dx(k,L)=-xl2+(-xl2-dx(k,L))
!		 elseif((yl2-0.02)<dy(k,L).and.dy(k,L)<yl2.and.(-xl2+0.02)<=dx(k,L).and.dx(k,L)<(xl2-0.02))then
!		   dy(k,L)=yl2-dy(k,L)+yl2
!		 elseif(-yl2<dy(k,L).and.dy(k,L)<(-yl2+0.02).and.(-xl2+0.02)<=dx(k,L).and.dx(k,L)<(xl2-0.02))then
!		   dy(k,L)=-yl2+(-yl2-dy(k,L))
!		 elseif((xl2-0.02)<dx(k,L).and.dx(k,L)<xl2.and.(yl2-0.02)<dy(k,L).and.dy(k,L)<yl2)then
!		   dy(k,L)=yl2-dy(k,L)+yl2
!		   dx(k,L)=xl2-dx(k,L)+xl2
!		 elseif((xl2-0.02)<dx(k,L).and.dx(k,L)<xl2.and.-yl2<dy(k,L).and.dy(k,L)<(-yl2+0.02))then
!		   dy(k,L)=-yl2+(-yl2-dy(k,L))
!		   dx(k,L)=xl2-dx(k,L)+xl2
!		 elseif(-xl2<dx(k,L).and.dx(k,L)<(-xl2+0.02).and.(yl2-0.02)<dy(k,L).and.dy(k,L)<yl2)then
!		   dy(k,L)=yl2-dy(k,L)+yl2
!		   dx(k,L)=-xl2+(-xl2-dx(k,L))
!		 elseif(-xl2<dx(k,L).and.dx(k,L)<(-xl2+0.02).and.-yl2<dy(k,L).and.dy(k,L)<(-yl2+0.02))then
!		   dx(k,L)=-xl2+(-xl2-dx(k,L))
!		   dy(k,L)=-yl2+(-yl2-dy(k,L))
!		 else
!		 continue
!		 end if
!	  end do
 ! end do

!--------------------受力分析-------------------
!----------原点涡坐标代换及生成新点涡的初始坐标---------
   do L=1,knp
      do k=1,kk
	     Gx(k+1,L)=dx(k,L)
		 Gy(k+1,L)=dy(k,L)
	  end do 
   end do
   do i=1,knp
      k=nt(i)
	  Gx(1,i)=bv(k,1)
	  Gy(1,i)=bv(k,2)
   end do

   !do m=1,mpv
      do i=1,knp
         self=nt(i)
	     spv(i,1)=0
	     spv(i,2)=0
	     do i2=1,nall
	        if(self==i2)then
		      continue
		    else
		      dis=(Gx(1,i)-bv(i2,1))**2+(Gy(1,i)-bv(i2,2))**2
		      spv(i,1)=spv(i,1)+btao(i2)*(-(Gy(1,i)-bv(i2,2)))/2/pi/dis
		      spv(i,2)=spv(i,2)+btao(i2)*(Gx(1,i)-bv(i2,1))/2/pi/dis
		    end if
	     end do
	     do L=1,knp
	        do k=1,kk
		       if(ndp(k,L)==0)then 
			     continue
		       else
			     dis=sqrt((Gx(1,i)-dx(k,L))**2+(Gy(1,i)-dy(k,L))**2)
			     if(dis<sigma(k))then
			       dis=sigma(k)
			     else
			       continue
			     end if
			     spv(i,1)=spv(i,1)+tao(k,L)*(-(Gy(1,i)-dy(k,L)))/2/pi/dis/dis
			     spv(i,2)=spv(i,2)+tao(k,L)*(Gx(1,i)-dx(k,L))/2/pi/dis/dis
			   end if
		    end do
	     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)*dt
	     Gy(1,i)=Gy(1,i)+spv(i,2)*dt
      end do
   end do
   do L=1,knp
      do k=1,kk
	  k1=kk+1-k
	  tao(k1+1,L)=tao(k1,L)
	  ndp(k1+1,L)=ndp(k1,L)
	  end do
   end do
   do i=1,knp
      k=nt(i)
	  ndp(1,i)=1
	  tao(1,i)=btao(k)
   !end do
 nstep=nstep+1

!end do
 open(unit=10,file='dx.txt')
 open(unit=20,file='dy.txt')
 write(10,"(1f8.4)")((dx(i,j),i=1,nlimit),j=1,knp)
 write(20,"(1f8.4)")((dy(i,j),i=1,nlimit),j=1,knp)
!write(*,"(8f8.4)")(btao(i),i=1,nall)
 close(10)
 close(20)
  !write(*,*)((ndp(i,j),i=1,nlimit),j=1,knp)