📄 dvm_circle.f90
字号:
subroutine circle()
use varible
logical:: state
real :: xl2,yl2,tkk,skk,ubb,tempv,dis
real :: bv0(nmax,2),bc0(nall,2),spv(knp,2)
integer:: i,k,j,L,M,k1,k2,L1,self,i2
!---------dx,dy为GX，GY数据存贮，FA为边界涡的数据存贮----
real :: dx(nlimit,knp),dy(nlimit,knp),VAx(nlimit,knp),VAy(nlimit,knp),FA(nall)
real :: bphy(nmax),kk1(6),kk2(6),kk3(6),kk4(6),fvx1(pnum),fvx2(pnum),fvy1(pnum),fvy2(pnum)
real ::px0(nlimit,pnum),py0(nlimit,pnum),pz0(nlimit,pnum),pvx0(nlimit,pnum),pvy0(nlimit,pnum),pvz0(nlimit,pnum)
!---------------------颗粒变量声明-----------------------
real ::bsigma,tvx,tvy,cosa,sina
!-----------参量赋值---------
open(unit=110,file="kk.txt")
xl2=xl/2
yl2=yl/2
ubb=0.5*U*U*yl
ddx=0.01
ddy=0.01
dnt=dt*U/xl
do k=kk+1,ns
fx(k)=0.0
fy(k)=0.0
end do
nstep=1
do i=1,nall
do j=1,2
bv0(i,j)=bv(i,j)
bc0(i,j)=bc(i,j)
end do
end do
!--------------------------------
state=.false.
!---------循环主体-------------
do while(kk<ns)
kk=kk+1
tkk=float(kk)*dt
if (state==.true.)then
skk=5.e-3*cos(2.0*pi*14.8*tkk)
do i=1,nall
bv(i,2)=bv0(i,2)+5.e-3*cos(2.0*pi*14.8*tkk)
bc(i,2)=bc0(i,2)+5.e-3*cos(2.0*pi*14.8*tkk)
end do
else
continue
end if
!--------------------------write-------------------------
!-------------KK步以前的涡的合并与坐标更新---------------
call combination(kk,Gx,Gy)
!-------------------已有数据存贮-------------------------
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)*(-sin(thta(i))*(bc(i,1)-dx(k,L))+cos(thta(i))*(bc(i,2)-dy(k,L)))/2.0/pi/dis/dis
end if
end do
end do
bphy(i)=-cos(thta(i))*U*cos(alpha)-sin(thta(i))*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)
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!------------计算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
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)=VAy(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)-VAy(k,L)
end if
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)=dy(k,L)+VAy(k,L)*ddt
end if
end do
end do
call combination(kk,dx,dy)
end do
!**********************************************************************
!**********************************************************************
!**********************************************************************
!与壁面接触时
do k=1,kk
do L=1,knp
dis=sqrt(dx(k,L)**2+dy(k,L)**2)
if(dis<xl2)then
tao(k,L)=0.0
dx(k,L)=0.0
dy(k,L)=0.0
ndp(k,L)=0.0
else
continue
end if
end do
end do
!----------------受力分析--------------
if (kk>=2)then
do i=1,nall
Fx(kk-1)=Fx(kk-1)+bv(i,2)*(btao(i)-Fa(i))/dt/ubb
Fy(kk-1)=Fy(kk-1)+bv(i,1)*(btao(i)-Fa(i))/dt/ubb
end do
else
continue
end if
!****************************************************************************
!****************************************************************************
!****************************************************************************
!****************************************************************************
if (kk==ns)then
continue
else
do k=1,kk
do L=1,knp
if (ndp(k,L)==0)then
continue
else
Fx(kk)=Fx(kk)+tao(k,L)*(dy(k,L)-Gy(k,L))/dt/Ubb
Fy(kk)=Fy(kk)+tao(k,L)*(dx(k,L)-Gx(k,L))/dt/ubb
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)*(1+pi/nall)
Gy(1,i)=bv(k,2)*(1+pi/nall)
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
do L=1,knp
Fx(kk)=Fx(kk)+Gy(1,L)*tao(1,L)/dt/ubb
Fy(kk)=Fy(kk)+Gx(1,L)*tao(1,L)/dt/ubb
end do
Fy(kk)=Fy(kk)+0.01*5.e-3*(2*pi*14.8)**2*cos(2*pi*14.8*tkk)/ubb
end if
!*************************************************************************************
!*************************************************************************************
!-----------------颗粒相的计算------------------------
write(110,*)kk
if(kk<=pt)then
continue
else
do i=1,pnum
fvx(2,i)=0.0
fvy(2,i)=0.0
do k=1,kk+1
do L=1,knp
if(ndp(k,L)==0)then
continue
else
dis=sqrt((xL*px(kk-1,i)-Gx(k,L))**2+(xL*py(kk-1,i)-Gy(k,L))**2)
if(dis<sigma(k))then
dis=sigma(k)
else
continue
end if
fvx(2,i)=fvx(2,i)+tao(k,L)*(-(xL*py(kk-1,i)-Gy(k,L)))/2/pi/dis/dis
fvy(2,i)=fvy(2,i)+tao(k,L)*(xL*px(kk-1,i)-Gx(k,L))/2/pi/dis/dis
end if
end do
end do
do j=1,nall
if(bnp(j)==1)then
continue
else
dis=sqrt((xL*px(kk-1,i)-bv(j,1))**2+(xL*py(kk-1,i)-bv(j,2))**2)
bsigma=max(abs(0.5*(bv(j+1,1)-bv(j,1))),abs(0.5*(bv(j+1,2)-bv(j,2))))
if(dis<bsigma)then
dis=bsigma
else
continue
end if
fvx(2,i)=fvx(2,i)+btao(j)*(-(xL*py(kk-1,i)-bv(j,2)))/2/pi/dis/dis
fvy(2,i)=fvy(2,i)+btao(j)*(xL*px(kk-1,i)-bv(j,1))/2/pi/dis/dis
end if
end do
fvx(2,i)=U*cos(alpha)-fvx(2,i)
fvy(2,i)=U*sin(alpha)-fvy(2,i)
!------------------------------------------------------------------------------------
dudx(i)=0.0
dudy(i)=0.0
dvdx(i)=0.0
dvdy(i)=0.0
do k=1,kk+1
do L=1,knp
if(ndp(k,L)==0)then
continue
else
dis=sqrt((xL*px(kk-1,i)-Gx(k,L))**2+(xL*py(kk-1,i)-Gy(k,L))**2)
if(dis<sigma(k))then
dis=sigma(k)
else
continue
end if
dudx(i)=dudx(i)-tao(k,L)*(xL*px(kk-1,i)-Gx(k,L))*(xL*py(kk-1,i)-Gy(k,L))/(dis**4)/pi
dudy(i)=dudy(i)+tao(k,L)*(dis**2-2*(xL*py(kk-1,i)-Gy(k,L))**2)/(dis**4)/pi/2.0
dvdx(i)=dvdx(i)-tao(k,L)*(dis**2-2*(xL*px(kk-1,i)-Gx(k,L))**2)/(dis**4)/pi/2.0
dvdy(i)=-dudx(i)
end if
end do
end do
fvx(2,i)=fvx(2,i)/U
fvy(2,i)=fvy(2,i)/U
dudx(i)=dudx(i)*Xl/U
dudy(i)=dudy(i)*XL/U
dvdx(i)=dvdx(i)*XL/U
dvdy(i)=dvdy(i)*XL/U
!-----------------------颗粒运动计算--------------------------
do j=1,6
kk1(j)=dfunc(ti,pvx(kk-1,i),pvy(kk-1,i),pvz(kk-1,i),px(kk-1,i),py(kk-1,i),pz(kk-1,i),j,i)
end do
do j=1,6
kk2(j)=dfunc(ti+dnt/2,pvx(kk-1,i)+kk1(1)*dnt/2,pvy(kk-1,i)+kk1(2)*dnt/2,pvz(kk-1,i)+kk1(3)*dnt/2,px(kk-1,i)+kk1(4)*dnt/2,py(kk-1,i)+kk1(5)*dnt/2,pz(kk-1,i)+kk1(6)*dnt/2,j,i)
end do
do j=1,6
kk3(j)=dfunc(ti+dnt/2,pvx(kk-1,i)+kk2(1)*dnt/2,pvy(kk-1,i)+kk2(2)*dnt/2,pvz(kk-1,i)+kk2(3)*dnt/2,px(kk-1,i)+kk2(4)*dnt/2,py(kk-1,i)+kk2(5)*dnt/2,pz(kk-1,i)+kk2(6)*dnt/2,j,i)
end do
do j=1,6
kk4(j)=dfunc(ti+dnt/2,pvx(kk-1,i)+kk3(1)*dnt/2,pvy(kk-1,i)+kk3(2)*dnt/2,pvz(kk-1,i)+kk3(3)*dnt/2,px(kk-1,i)+kk3(4)*dnt/2,py(kk-1,i)+kk3(5)*dnt/2,pz(kk-1,i)+kk3(6)*dnt/2,j,i)
end do
pvx0(kk,i)=pvx(kk-1,i)
pvy0(kk,i)=pvy(kk-1,i)
pvz0(kk,i)=pvz(kk-1,i)
px0(kk,i)=px(kk-1,i)
py0(kk,i)=py(kk-1,i)
pz0(kk,i)=pz(kk-1,i)
pvx(kk,i)=pvx(kk-1,i)+dnt*(kk1(1)+2*kk2(1)+2*kk3(1)+kk4(1))/6
pvy(kk,i)=pvy(kk-1,i)+dnt*(kk1(2)+2*kk2(2)+2*kk3(2)+kk4(2))/6
pvz(kk,i)=pvz(kk-1,i)+dnt*(kk1(3)+2*kk2(3)+2*kk3(3)+kk4(3))/6
px(kk,i)=px(kk-1,i)+dnt*(kk1(4)+2*kk2(4)+2*kk3(4)+kk4(4))/6
py(kk,i)=py(kk-1,i)+dnt*(kk1(5)+2*kk2(5)+2*kk3(5)+kk4(5))/6
pz(kk,i)=pz(kk-1,i)+dnt*(kk1(6)+2*kk2(6)+2*kk3(6)+kk4(6))/6
!=====================================================================
fvx(1,i)=fvx(2,i)
fvy(1,i)=fvy(2,i)
!=====================================================================
dis=xL*sqrt(px(kk,i)**2+py(kk,i)**2)
dis0=xl2/dis
dis1=xl2/sqrt(px0(kk,i)**2+py0(kk,i)**2)/xL
if(dis>xL2)then
continue
else
cosa=xl*px(kk,i)/dis
sina=xl*py(kk,i)/dis
Px(kk,i)=0.5*(px0(kk,i)*dis1+px(kk,i)*dis0)
Py(kk,i)=0.5*(py0(kk,i)*dis1+px(kk,i)*dis0)
tvx=pvx(kk,i)
tvy=pvy(kk,i)
pvx(kk,i)=-(tvx*cosa+tvy*sina)*cosa+(tvx*sina-tvy*cosa)*sina
pvy(kk,i)=-(tvx*cosa+tvy*sina)*sina-(tvx*sina-tvy*cosa)*cosa
end if
end do
ti=ti+dnt
!-----------------颗粒边界条件-------------------
end if
write(*,*)"kk=",kk,"ti=",ti
!write(*,*)(dvdy(i),i=1,pnum)
end do
DO I=1,NLIMIT
DO J=1,PNUM
PVX(I,J)=PVX(I,J)*U
PVY(I,J)=PVY(I,J)*U
PX(I,J)=PX(I,J)*XL
PX(I,J)=PX(I,J)*XL
END DO
END DO
close(110)
end subroutine
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -