[Fortran] 纯文本查看 复制代码
program DVM
!
implicit none !-->all the variables MUST be declared
!
integer,parameter :: n = 10
!Declaration of variable
real(8) :: U,a,dt
real(8) :: xs,ys
integer :: p,q
integer :: i
!p for which vortex it is, q for how long it exist
real(8),dimension(n-1):: gama
!the circulation of each vortex
real(8),dimension(n,n) :: xx,yy,xxi,yyi
!the location of the vortices and their corresponding images, the first index refers to which vortex it is, and the second one refers to how long it have existed
print*,'input the Velocity of the free stream, the Radius of the cylinder'
read*,U,a
print*,'input the Unit Time Step'
read*,dt
xs = 0
ys = a
!uu = U(1-((xx(p,q)**2 - yy(p,q)**2)*a**2)/((xx(p,q)**2 + yy(p,q)**2)**2))
!vv = -2*U*xx(p,q)*yy(p,q)*a**2/(xx(p,q)**2+yy(p,q)**2)**2
!formula (14) (15) in the lecture note
do p = 1,n-1
!set q = 1 to get the new pop-up vortex---------------------
i = 1
q = 1
!define the value of x of every pop-up vortex---------------
xx(p,q) = 0
!define the value of y of every pop-up vortex--------------
if (MOD(p,2)/=0) then
yy(p,q) = 1.1*a
else
yy(p,q) = -1.1*a
end if
!calculate gamma of the pop-up vortex----------------------
q = p
call getgama()
q = p
!calculate the new location of the vortices at present-----
do i = 1,p
call getlocation()
q = q - 1
end do
end do
print*,gama
contains
subroutine image(xxx,yyy,xxxi,yyyi)
implicit none
real(8) :: xxx,yyy,xxxi,yyyi
xxxi = a**2*(xxx/(xxx**2 + yyy**2))
yyyi = a**2*(yyy/(xxx**2 + yyy**2))
return
end subroutine image
subroutine getgama()
implicit none
real(8), dimension(p) :: simp
!to simplify the calculation, use simp to represent the iterative part of formula (18)
real(8), parameter :: Pi = 3.1415927
if(p==1) then
simp(1) = -2*Pi*U*(1-((xs**2 - ys**2)*(a**2))/((xs**2 + ys**2)**2))
call image(xx(1,1),yy(1,1),xxi(1,1),yyi(1,1))
gama(1) = simp(1)/(-(ys-yy(1,1))/((xs-xx(1,1))**2+(ys-yy(1,1))**2)+(ys-yyi(1,1))/((xs-xxi(1,1))**2+(ys-yyi(1,1))**2))
else
if (MOD(p,2)/=0) then
do while (i<p)
call image(xx(i,q),yy(i,q),xxi(i,q),yyi(i,q))
simp(i+1) = simp(i)+gama(i)*(ys-yy(i,q))/((xs-xx(i,q))**2+(ys-yy(i,q))**2)-gama(i)*(ys-yyi(i,q))/((xs-xxi(i,q))**2+(ys-yyi(i,q))**2)
q = q-1
i = i+1
end do
gama(p) = simp(p)/(-(ys-yy(p,1))/((xs-xx(p,1))**2+(ys-yy(p,1))**2)+(ys-yyi(p,1))/((xs-xxi(p,1))**2+(ys-yyi(p,1))**2))
else
do while (i<p)
call image(xx(i,q),yy(i,q),xxi(i,q),yyi(i,q))
simp(i+1) = simp(i)+gama(i)*(-ys-yy(i,q))/((xs-xx(i,q))**2+(-ys-yy(i,q))**2)-gama(i)*(-ys-yyi(i,q))/((xs-xxi(i,q))**2+(-ys-yyi(i,q))**2)
q = q-1
i = i+1
end do
gama(p) = simp(p)/(-(-ys-yy(p,1))/((xs-xx(p,1))**2+(-ys-yy(p,1))**2)+(-ys-yyi(p,1))/((xs-xxi(p,1))**2+(-ys-yyi(p,1))**2))
end if
end if
return
end subroutine getgama
subroutine getlocation()
implicit none
!p tells how many vortices do we have now
!we should move (p-1) vortices
!q tells which stage is each vortex at
real, parameter :: Pi = 3.1415927
real(8) :: Ustream,Vstream
real(8) :: simpxm = 0,simpxp = 0,simpym = 0,simpyp = 0
integer :: j,m
m = q
do j = 1,p
if(j/=i) then
q = p + 1 - j
simpxm = simpxm + dt*( Gama(j) / (2*Pi) ) * (yy(i,q) - yy(j,q)) / ((xx(i,q)-xx(j,q))**2 + (yy(i,q)-yy(j,q))**2)
end if
end do
do j = 1,p
q = p + 1 - j
simpxp = simpxp + dt*( Gama(j) / (2*Pi) ) * (yy(i,q) - yyi(j,q)) / ((xx(i,q)-xxi(j,q))**2 +(yy(i,q)-yyi(j,q))**2)
end do
do j = 1,p
if(j/=i) then
q = p + 1 - j
simpym = simpym + dt*( Gama(j) / (2*Pi) ) * (xx(j,q) - xx(i,q)) / ((xx(i,q)-xx(j,q))**2 + (yy(i,q)-yy(j,q))**2)
end if
end do
do j = 1,p
q = p + 1 - j
simpyp = simpyp + dt*( Gama(j) / (2*Pi) ) * (xxi(j,q) - xx(i,q)) / ((xx(i,q)-xxi(j,q))**2 +(yy(i,q)-yyi(j,q))**2)
end do
Ustream = U * (1-(xx(i,m)**2-yy(i,m)**2) * a**2 / (xx(i,m)**2+yy(i,m)**2)**2)
Vstream = -2 * U * xx(i,m) * yy(i,m) * a**2 / (xx(i,m)**2+yy(i,m)**2)**2
xx(i,m+1) = xx(i,m) + Ustream * dt - simpxm + simpxp
yy(i,m+1) = yy(i,m) + Vstream * dt - simpym + simpyp
return
end subroutine getlocation
end program DVM
这是留的大作业，研究好久了，感觉不难，但做起来很繁琐，麻烦大家了，另外请问有没有调试程序的好方法，一个一个找太累了！！！