dvm_setup.f90

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

subroutine setup()
use varible
real::deltx,delty,dis
real::bphy(nmax)
integer i,j,k,L,K1,L1
!--------------物理参数设置----------------------
dudxl=U/xl
alpha=0.0
alpha=alpha*pi/180.0
dt=0.001
re=U*xl/fmiu
ddt=dt/float(Mpv)
Lc=pi*xl          
dL=Lc/nall              !GD
dthta=2*pi/nall
sigma1=dL/4.0
!------------------------------------------------
KK=1
nstep=1
ddt=dt/float(mpv)
call nsp()
!------------几何参数设置------------------------
do L=1,knp
   do k=1,Ns+1
   Ndp(k,L)=1
   end do
end do
!------------壁面点涡坐标设置--------------------
do i=1,nall
bv(i,1)=xl/2*cos((i-1)*dthta)
bv(i,2)=xl/2*sin((i-1)*dthta)
thta(i)=(i-1)*dthta
end do
bv(nmax,1)=bv(1,1)
bv(nmax,2)=bv(1,2)

!--------------壁面控制点坐标设置----------------
do i=1,nall
bc(i,1)=(bv(i,1)+bv(i+1,1))/2.0
bc(i,2)=(bv(i,2)+bv(i+1,2))/2.0
end do 
                                                                                                                                                                                                                                                                                                                                                                                                                                                                   
!-------------壁面计算系数设置-------------------



!------------------dphy(i)-----------------------
do i=1,nall
bphy(i)=-cos(thta(i))*U*cos(alpha)-sin(thta(i))*U*sin(alpha)
end do
bphy(nmax)=0.0

!------------------bmatrx(i,j)---------------------
do i=1,nall
  do j=1,nall
  dis=(bc(i,1)-bv(j,1))**2+(bc(i,2)-bv(j,2))**2
  matrx1(i,j)=1/2.0/pi/dis*(cos(thta(i))*(bc(i,2)-bv(j,2))-sin(thta(i))*(bc(i,1)-bv(j,1)))
  end do
end do
!-----------------去掉一个多余条件-----------------
do j=1,nall
matrx1(nmax,j)=1.0
end do
bphy(nmax)=0.0
do i=1,nall
   do j=1,nall
   bmatrx(i,j)=0.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)

open(unit=80,file='bv10.txt')
open(unit=90,file='bv20.txt')
write(80,"(1f8.4)")(bv(i,1),i=1,nall)
write(90,"(1f8.4)")(bv(i,2),i=1,nall)

!******************************************************
!-------------颗粒初始条件设置-------------------------
ti=0
pt=20
do i=1,pnum
   px(pt,i)=-0.08/xl
   py(pt,i)=(-0.1+(i-1)*0.004)/xl
   pvx(pt,i)=5.0/U
   pvy(pt,i)=0.0/U
   pz(pt,i)=1.0
   pvz(pt,i)=0.1
end do 

do j=1,pnum
   fvx(1,j)=0.0
   fvy(1,j)=0.0
end do
return
end subroutine setup