📄 youngs-vof.f90
字号:
! calculate the transportation from neighbour cell
subroutine transport(c,itype, tanalfa, tanbeta, u1,u3,u4,u2,dx,dy,dt,f1,f3,f4,f2)
implicit double precision (a-h,o-z)
cotalfa=1.0/tanalfa
cotbeta=1.0/tanbeta
f1=0
f3=0
f4=0
f2=0
if(itype.eq.1)then
s1=0
s3=sqrt(2.0*c*tanalfa)
s4=0
s2=sqrt(2.0*c*tanalfa)
if(u1.gt.0)then
if(u1*dt.le.(1.0-s2)*dy)then
f1=0
else
temp1=u1*dt-(1.0-s2)*dy
f1=0.5*temp1*temp1*cotbeta
endif
endif
if(u2.gt.0)then
if(u2*dt.le.s3*dx)then
f2=c*dx*dy
else
f2=0.5*u2*dt*(2.0-u2*dt/(s3*dx))*s2*dy
endif
endif
if(u3.lt.0)then
if((abs(u3)*dt.ge.s2*dy))then
f3=c*dx*dy
else
f3=0.5*abs(u3)*dt*(2.0-abs(u3)*dt/(s2*dy))*s3*dx
endif
endif
if(u4.lt.0)then
if(abs(u4)*dt.le.(1.0-s3)*dx)then
f4=0
else
temp1=abs(u4)*dt-(1.0-s3)*dx
f4=0.5*temp1*temp1*tanbeta
endif
endif
else if(itype.eq.2)then
s1=0
s3=1.0
s4=c-0.5*tanalfa
s2=c+0.5*tanalfa
if(u1.gt.0)then
if(u1*dt.le.(1.0-s2)*dy)then
f1=0
else if(u1*dt.le.(1.0-s4)*dy)then
temp1=u1*dt-(1.0-s2)*dy
f1=0.5*temp1*temp1*cotbeta
else
f1=u1*dt*dx-(1.0-c)*dx*dy
endif
endif
if(u2.gt.0)then
f2=u2*dt*(s2*dy-0.5*u2*dt*tanbeta)
endif
if(u3.lt.0)then
if((abs(u3)*dt.le.s4*dy))then
f3=abs(u3)*dt*dx
else if((abs(u3)*dt.le.s2*dy))then
temp1=abs(u3)*dt-s4*dy
f3=abs(u3)*dt*dx-0.5*temp1*temp1*cotbeta
else
f3=c*dx*dy
endif
endif
if(u4.lt.0)then
f4=abs(u4)*dt*(s4*dy+0.5*abs(u4)*dt*tanbeta)
endif
else if(itype.eq.3)then
s1=c-0.5*tanalfa
s3=c+0.5*tanalfa
s4=0
s2=1.0
if(u1.gt.0)then
f1=u1*dt*(s1*dx+0.5*u1*dt*cotbeta)
endif
if(u2.gt.0)then
if(u2*dt.le.s1*dx)then
f2=u2*dt*dy
else if(u2*dt.le.s3*dx)then
temp1=u2*dt-s1*dx
f2=u2*dt*dy-0.5*temp1*temp1*tanbeta
else
f2=c*dx*dy
endif
endif
if(u3.lt.0)then
f3=abs(u3)*dt*(s3*dy+0.5*abs(u3)*dt*cotbeta)
endif
if(u4.lt.0)then
if((abs(u4)*dt.le.(1.0-s3)*dx))then
f4=0
else if((abs(u4)*dt.le.(1.0-s1)*dx))then
temp1=abs(u4)*dt-(1.0-s3)*dx
f4=0.5*temp1*temp1*tanbeta
else
f4=abs(u4)*dt*dy-(1.0-c)*dx*dy
endif
endif
else if(itype.eq.4)then
s1=1.0-sqrt(2.0*(1.0-c)*cotalfa)
s3=1.0
s4=1.0-sqrt(2.0*(1.0-c)*tanalfa)
s2=1.0
if(u1.gt.0)then
if(u1*dt.ge.(1.0-s4)*dy)then
f1=u1*dt*dx-(1.0-c)*dx*dy
else
f1=u1*dt*(s1*dx+0.5*u1*dt*cotbeta)
endif
endif
if(u2.gt.0)then
if(u2*dt.le.s1*dx)then
f2=u2*dt*dy
else
temp1=u2*dt-s1*dx
f2=u2*dt*dy-0.5*temp1*temp1*tanbeta
endif
endif
if(u3.lt.0)then
if(abs(u3)*dt.le.s4*dy)then
f3=abs(u3)*dt*dx
else
temp1=abs(u3)*dt-s4*dy
f3=abs(u3)*dt*dx-0.5*temp1*temp1*cotbeta
endif
endif
if(u4.lt.0)then
if((abs(u4)*dt.ge.(1.0-s1)*dx))then
f4=abs(u4)*dt*dy-(1.0-c)*dx*dy
else
f4=abs(u4)*dt*(s4*dy+0.5*abs(u4)*dt*tanbeta)
endif
endif
else
write(*,*) 'error in subroutine transport'
stop
endif
end
! program main
! solve fluid volume function equation using Young's VOF method
program main
parameter(isize=100,pi=3.14159265)
implicit double precision (a-h,o-z)
double precision c0(0:isize,0:isize),c1(0:isize,0:isize)
double precision u(0:isize,0:isize),v(0:isize,0:isize)
double precision x(0:isize),y(0:isize)
open(unit=1,file='d:/you1_0a.plt',status='unknown')
open(unit=2,file='d:/you1_0b.plt',status='unknown')
write(*,*) 'Input T:'
read(*,*) supt
h=1.0/isize
dt=0.1*h
dx=h
dy=h
emikro=1.0e-10
emk2=0
emk2=0
x0=0.5
y0=0.3
r1=0.2
r2=0.15
do i=0,isize
x(i)=dx*i
y(i)=dy*i
enddo
do i=0,isize
do j=0,isize
u(i,j)=-pi*cos(pi*(x(i)-0.5))*sin(pi*(y(j)-0.5))
v(i,j)=pi*sin(pi*(x(i)-0.5))*cos(pi*(y(j)-0.5))
enddo
enddo
do i=1,isize
do j=0,isize
! if(dist((x(i),y(j)),(x0,y0)).le.r1)then
if(sqrt((x(i)-x0)*(x(i)-x0)+(y(j)-y0)*(y(j)-y0)).le.r1)then
c0(i,j)=1.0
else
c0(i,j)=0
endif
enddo
enddo
do ll=1,2
t=0
it=0
do while(t.lt.supt)
do i=1,isize
do j=1,isize
c1(i,j)=c0(i,j)
enddo
enddo
do i=1,isize-1
do j=1,isize-1
ft=0
fb=0
fl=0
fr=0
c=c0(i,j)
ut=v(i,j)
ub=v(i,j-1)
ul=u(i-1,j)
ur=u(i,j)
if (c.ge.1.0-emikro)then
if(ut.gt.0)then
ft=ut*dt*dx*c
endif
if(ub.lt.0)then
fb=abs(ub)*dt*dx*c
endif
if(ul.lt.0)then
fl=abs(ul)*dt*dy*c
endif
if(ur.gt.0)then
fr=ur*dt*dy*c
endif
else if(c.gt.0)then
rnx=(c0(i+1,j+1)+2.0*c0(i+1,j)+c0(i+1,j-1)-c0(i-1,j+1)-2.0*c0(i-1,j)-c0(i-1,j-1))/dx
rny=(c0(i+1,j+1)+2.0*c0(i,j+1)+c0(i-1,j+1)-c0(i+1,j-1)-2.0*c0(i,j-1)-c0(i-1,j-1))/dy
if(abs(rny).le.emk2)then
if(rnx.ge.0)then
u1=ut
u2=ur
u3=ub
u4=ul
else
u1=ut
u2=-ul
u3=ub
u4=-ur
endif
f1=0
f2=0
f3=0
f4=0
if(u1.gt.0)then
f1=abs(u1)*dt*c*dx
endif
if(u3.lt.0)then
f3=abs(u3)*dt*c*dx
endif
if(u2.gt.0)then
if(abs(u2)*dt.le.c*dx)then
f2=abs(u2)*dt*dy
else
f2=c*dx*dy
endif
endif
if(u4.lt.0)then
if(abs(u4)*dt.le.(1.0-c)*dx)then
f4=0
else
f4=(abs(u4)*dt-(1.0-c)*dx)*dy
endif
endif
if(rnx.ge.0)then
ft=f1
fb=f3
fr=f2
fl=f4
else
ft=f1
fb=f3
fr=f4
fl=f2
endif
else if(abs(rnx).le.emk2)then
if(rny.ge.0)then
u1=ut
u2=ur
u3=ub
u4=ul
else
u1=-ub
u2=ur
u3=-ut
u4=ul
endif
f1=0
f2=0
f3=0
f4=0
if(u1.gt.0)then
if(u1*dt.le.c*dy)then
f1=u1*dt*dx
else
f1=c*dx*dy
endif
endif
if(u3.lt.0)then
if(abs(u3)*dt.le.(1.0-c)*dy)then
f3=0
else
f3=(abs(u3)*dt-(1.0-c)*dy)*dx
endif
endif
if(u2.gt.0)then
f2=u2*dt*c*dy
endif
if(u4.lt.0)then
f4=abs(u4)*dt*c*dy
endif
if(rnx.ge.0)then
ft=f1
fb=f3
fr=f2
fl=f4
else
ft=f3
fb=f1
fr=f2
fl=f4
endif
else
if(rnx.gt.0.and.rny.lt.0)then
u1=ut
u2=ur
u3=ub
u4=ul
else if(rnx.lt.0.and.rny.gt.0)then
u1=-ub
u2=-ul
u3=-ut
u4=-ur
else if(rnx.lt.0.and.rny.lt.0)then
u1=ut
u2=-ul
u3=ub
u4=-ur
else
u1=-ub
u2=ur
u3=-ut
u4=ul
endif
rnx1=abs(rnx)
rny1=-abs(rny)
tanbeta=-rnx1/rny1
cotbeta=1.0/tanbeta
tanalfa=dx/dy*tanbeta
cotalfa=1.0/tanalfa
if(tanalfa.le.1.0)then
if(c.le.0.5*tanalfa)then
itype=1
else if(c.le.1.0-0.5*tanalfa)then
itype=2
else
itype=4
endif
else
if(c.le.0.5*cotalfa)then
itype=1
else if(c.le.1.0-0.5*cotalfa)then
itype=3
else
itype=4
endif
endif
call transport(c,itype, tanalfa, tanbeta, u1,u3,u4,u2,dx,dy,dt,f1,f3,f4,f2)
if(rnx.gt.0.and.rny.lt.0)then
ft=f1
fb=f3
fl=f4
fr=f2
else if(rnx.lt.0.and.rny.gt.0)then
ft=f3
fb=f1
fl=f2
fr=f4
else if(rnx.lt.0.and.rny.lt.0)then
ft=f1
fb=f3
fl=f2
fr=f4
else
ft=f3
fb=f1
fl=f4
fr=f2
endif
endif
endif
c1(i,j+1)=c1(i,j+1)+ft/dx/dy
c1(i,j-1)=c1(i,j-1)+fb/dx/dy
c1(i+1,j)=c1(i+1,j)+fr/dx/dy
c1(i,j)=c1(i,j)-(ft+fb+fl+fr)/dx/dy
enddo
enddo
do i=1,isize
do j=1,isize
c0(i,j)=max(min(c1(i,j),1.0),0.0)
enddo
enddo
t=t+1
it=it+1
enddo
write(ll,*) 'TITLE="EXANPLE:3D GEOMETRIES"'
write(ll,*) 'VARIABLES="X","Y","Z"'
write(ll,*) 'ZONE T="Floor",I=',isize,'J=',isize,'F=POINT'
do i=1,isize
do j=1,isize
write(ll,102) x(i),y(j),c0(i,j)
enddo
enddo
do i=1,isize
do j=1,isize
u(i,j)=pi*cos(pi*(x(i)-0.5))*sin(pi*(y(j)-0.5))
v(i,j)=-pi*sin(pi*(x(i)-0.5))*cos(pi*(y(j)-0.5))
enddo
enddo
enddo
102 format(1x,e20.10,e20.10,e20.10)
close(1)
close(2)
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -