particles_rkfunction.f90

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

real function dfunc(x1,x2,x3,x4,x5,x6,x7,k,i)
!===x1-t,x2-vpx,x3-vpy,x4-m,x5-xp,x6-yp,x7-?,k-记数器===
use varible
real   ::x1,x2,x3,x4,x5,x6,x7
real   ::A,B,ttao,mf,mp,rouf,roup,ra,rp,wp
integer::k,i
real   ::gr
!---------颗粒计算参数-------------
rd=0.001
rouf=1000
roup=600
g=9.8
mp=4/3*pi*rd**3*roup
mf=4/3*pi*rd**3*rouf
Rp=rouf/(roup+0.5*rouf)
A=6*pi*rd*fmiu*xl/U/(mp+0.5*mf)
wp=(mp-mf)*g/6/pi/rd/fmiu/U
select case(k)
case(1)
!d(pvx)/dt
dfunc=A*(fvx(2,i)-x2)+((fvx(2,i)+0.5*x2)*dudx(i)+(fvy(2,i)+0.5*x3)*dudy(i))*Rp+1.5*Rp*(fvx(2,i)-fvx(1,i))/dnt
!dfunc=10*(x2-x3)

case(2)
!d(pvy)/dt
dfunc=A*(fvy(2,i)-x3+Wp)+((fvx(2,i)+0.5*x2)*dvdx(i)+(fvy(2,i)+0.5*x3)*dvdy(i))*Rp+1.5*Rp*(fvy(2,i)-fvy(1,i))/dnt
!dfunc=x2*x4+25*x2-x3

case(3)
!d(m)/dt
dfunc=0.0
!dfunc=-8/3*x4+x2*x3

case(4)
!d(px)/dt
dfunc=x2


case(5)
!d(py)/dt
dfunc=x3


case(6)
!d(?)/dt


dfunc=0
case default
write(*,*) "data error!"
end select
return
end function dfunc