Fortran 언어 - 자유 인터페이스 프로그램
9009 단어 수치 시뮬레이션LBM자유 인터페이스freesurface
PROGRAM dambreak
USE IFPORT
implicit none
! //LBM model
Real*8, parameter:: w(0:8) = (/4.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0/)
Real*8, parameter:: rhoA = 1.d0
Integer, parameter:: e(0:8,0:1)=(/(/0,1,0,-1, 0,1,-1,-1, 1/),(/0,0,1, 0,-1,1, 1,-1,-1/)/)
Integer, parameter:: inv(0:8) = (/0,3,4,1,2,7,8,5,6/)
Integer, parameter:: xDim=5900,yDim=920,xFluid=3400,yFluid=600
integer, parameter:: itslip(0:8) = (/0,1,2,3,2,5,6,6,5/)
integer, parameter:: ibslip(0:8) = (/0,1,4,3,4,8,7,7,8/)
integer, parameter:: irslip(0:8) = (/0,1,2,1,4,5,5,8,8/)
integer, parameter:: ilslip(0:8) = (/0,3,2,3,4,6,6,7,7/)
! Real*8 :: w(0:8),e(0:8,0:1),inv(0:8) !Real*8
Real*8 :: rho(1:xDim,1:yDim),mass(1:xDim,1:yDim),u(1:xDim,1:yDim,0:1)
Real*8 :: fEq(1:xDim,1:yDim,0:8),f(1:xDim,1:yDim,0:8),fpost(1:xDim,1:yDim,0:8)
Real*8 :: uxy,uSqr,F_eq,mass1,mass2,P_lt2ph
Real*8 :: afb,aft,afr,afl,adb,adt,adr,adl
Real*8 :: dx,dt,vis_phy,gra_phy,den_phy
Real*8 :: gravn,vis_lb,tau,omega,times,p_lb
Integer:: iflag(0:xDim+1,0:yDim+1)
Integer:: i,tStep,x,y,tMax,BCTYPE
Character ( len = 100 ) :: command_file_name = 'file_commands.txt'
integer::csys
! \\Deleting the former folder
csys = SYSTEM('C:\Windows\System32\cmd.exe /E:ON /V:ON /K "phan.bat" ')
!
CALL Readparameters(tMax,BCTYPE,dx,dt,vis_phy,gra_phy,den_phy)
print*,tMax,BCTYPE,dx,dt,vis_phy,gra_phy,den_phy
if (BCTYPE.eq.1) then !Non-slip BC
adr = 1.0d0; afr=1.d0 - adr
adb = 1.0d0; afb=1.d0 - adb
adt = 1.0d0; aft=1.d0 - adt
adl = 1.0d0; afl=1.d0 - adl
elseif(BCTYPE.eq.2) then !Slip BC
adr = 0.0d0; afr=1.d0 - adr
adb = 0.0d0; afb=1.d0 - adb
adt = 0.0d0; aft=1.d0 - adt
adl = 0.0d0; afl=1.d0 - adl
elseif(BCTYPE.eq.3) then !Hybrid BC
adr = 0.0d0; afr=1.d0 - adr
adb = 0.0d0; afb=1.d0 - adb
adt = 1.0d0; aft=1.d0 - adt
adl = 0.0d0; afl=1.d0 - adl
endif
! //Unit conversion
P_lt2ph=dx*dx/(dt*dt)
! print*,P_lt2ph
! pause
! stop
! //Calculating relaxation time
gravn=gra_phy*dt*dt/dx !
vis_lb=vis_phy*dt/(dx*dx) !viscosity
tau=3.d0*vis_lb+0.5d0
omega=1.995d0
! omega=1.d0/tau
print*,omega,tau
!
CALL iniCondi(rho,mass,u,iflag,f,fpost,xDim,yDim,xFluid,yFluid,w)
!
! goto 301
open (11, file='pressure.txt',status='unknown')
open (12, file='masstotal.txt',status='unknown')
! --------------------------------------------------------------------------------------------------
DO 300 tStep = 1, tMax
CALL massUpdate(fpost,rho,mass,iflag,xDim,yDim,e,inv)
CALL stream(f,fpost,mass,iflag,u,afb,aft,afr,afl,adb,adt,adr,adl,xDIm,yDim,rhoA,e,w,inv,ibslip,itslip,irslip,ilslip)
CALL collide(f,fpost,iflag,rho,u,xDim,yDim,omega,gravn,w,e,p_lb,mass,P_lt2ph)
CALL cellsUpdate(f,rho,u,iflag,mass,xDim,yDim,w,e)
! //Print out result
IF(mod(tStep,100) .eq. 0)then
print*,'Time Step = ',tStep
! CALL output(rho,u,iflag,mass,tStep,xDim, yDim)
ENDIF
times=dt*dble(tstep) ! dble , 。
write(11,*) times,p_lb
!
do x=1,xDim
do y=1,yDim
do i=0,8
fpost(x,y,i)=f(x,y,i) ! f fpost
enddo
enddo
enddo
! //Checking mass conservation
mass1=0.d0
mass2=0.d0
do x=1,xDim
do y=1,yDim
mass1=mass1+mass(x,y) !
if (iflag(x,y).eq.1.or.iflag(x,y).eq.2) mass2=mass2+mass(x,y)
enddo
enddo
if (isnan(mass1)) then
Print*,'THE NaN ERROR OCCURS'
PAUSE
STOP
endif
Write(12,116)tStep, mass1,mass2,mass1-mass2
116 format(i15,3f20.5)
!
300 CONTINUE
close(11)
!
301 continue
! call run_gnuplot ( command_file_name )
!
END PROGRAM dambreak
!-----------Equilibrium Function-------------------------------------------------------
Real*8 function F_eq(wi,rhoxy,uxy,uSqr)
implicit none
Real*8:: uxy,uSqr,wi,rhoxy
F_eq = wi*rhoxy*(1.d0+3.0d0*uxy+4.5d0*uxy*uxy-1.5d0*uSqr)
end function
!-------------------------------------------------------------------------------------
!-----Reading parameters--------------------------------------------------------------
SUBROUTINE Readparameters(tMax,BCTYPE,dx,dt,vis_phy,gra_phy,den_phy)
!
implicit none
Integer::tMax,BCTYPE
Real*8:: dx,dt,vis_phy,gra_phy,den_phy
!
open(1,file='PARAMETERS.txt')
read(1,*)
read(1,*)
read(1,*)tMax
read(1,*)
read(1,*)dx,dt
read(1,*)
read(1,*)BCTYPE
read(1,*)
read(1,*)vis_phy
read(1,*)
read(1,*)gra_phy
read(1,*)
read(1,*)den_phy
close(1)
END SUBROUTINE
!--------------------------------------------------------------------------------------
!------------------SUBROUTINE: iniCondi------------------------------------------------
! Function: Generating Initial conditions
!
SUBROUTINE iniCondi(rho,mass,u,iflag,f,fpost,xDim,yDim,xFluid,yFluid,w)
implicit none
!
Integer:: xDim,yDim,xFluid,yFluid,x,y,i
Integer:: iflag(0:xDim+1,0:yDim+1)
Real*8 :: rho(1:xDim,1:yDim),mass(1:xDim,1:yDim),u(1:xDim,1:yDim,0:1)
Real*8 :: f(1:xDim,1:yDim,0:8),fpost(1:xDim,1:yDim,0:8)
Real*8 :: w(0:8)
! //For fluid cells
do x=1,xFluid
do y=1,yFluid
rho(x,y) = 1.d0
iflag(x,y)= 1
mass(x,y) = rho(x,y)
u(x,y,0) = 0.d0
u(x,y,1) = 0.d0
do i = 0, 8
f(x,y,i)=w(i)
fpost(x,y,i)=w(i)
end do
enddo
enddo
! //For surface cells
y=yFluid+1 ! ， mass= ，
do x=1,xFluid
rho(x,y) = 1.d0
iflag(x,y)= 2
mass(x,y) = rho(x,y)*0.5d0
u(x,y,0) = 0.d0
u(x,y,1) = 0.d0
do i = 0, 8
f(x,y,i)=w(i)
fpost(x,y,i)=w(i)
end do
enddo
x=xFluid+1 ! ， mass=
do y=1,yFluid
rho(x,y) = 1.d0
iflag(x,y)= 2
mass(x,y) = rho(x,y)*0.5d0
u(x,y,0) = 0.d0
u(x,y,1) = 0.d0
do i = 0, 8
f(x,y,i)=w(i)
fpost(x,y,i)=w(i)
end do
enddo
x=xFluid+1! ，mass=
y=yFluid+1
rho(x,y) = 1.d0
iflag(x,y)= 2
mass(x,y) = rho(x,y)*0.25d0 !mass rho*(1+kapa) ！
u(x,y,0) = 0.d0
u(x,y,1) = 0.d0
do i = 0, 8
f(x,y,i)=w(i)
fpost(x,y,i)=w(i)
end do
! //For gas cells
do x=xFluid+2,xDim
do y=1,yDim
rho(x,y) = 0.d0
iflag(x,y)= 3
mass(x,y) = 0.d0
u(x,y,0) = 0.d0
u(x,y,1) = 0.d0
enddo
enddo
do x=1,xFluid+1
do y=yFluid+2,yDim
rho(x,y) = 0.d0
iflag(x,y)= 3
mass(x,y) = 0.d0
u(x,y,0) = 0.d0
u(x,y,1) = 0.d0
enddo
enddo
! //For wall boundary flags
do x=0,xDim+1,xDim+1 ! =xDim+1，
do y=0,yDim+1
iflag(x,y)= 0
enddo
enddo
do x=0,xDim+1
do y=0,yDim+1,yDim+1
iflag(x,y)= 0
enddo
enddo
END SUBROUTINE iniCondi
!--------------------------------------------------------------------------------------
!-----------------SUBROUTINE: massUpdate-----------------------------------------------
! Function: Evaluate mass-flux
!
SUBROUTINE massUpdate(fpost,rho,mass,iflag,xDim,yDim,e,inv)
implicit none
!
Integer:: xDim,yDim
Real*8 :: fpost(1:xDim,1:yDim,0:8),rho(1:xDim,1:yDim)
Real*8 :: mass(1:xDim,1:yDim)
Real*8 :: eps,epst,dmtol,dm
Integer:: e(0:8,0:1),inv(0:8)
Integer:: iflag(0:xDim+1,0:yDim+1)
Integer:: x,y,i,xt,yt
!
DO 270 x = 1, xDim
DO 270 y = 1, yDim
if (iflag(x,y).eq.1) then !For fluid
mass(x,y)=rho(x,y)
elseif (iflag(x,y).eq.2) then !For interface
dmtol=0.d0 !mass
eps = 0.d0 !eps 。
if (rho(x,y).gt.0.d0) eps = mass(x,y)/rho(x,y) !GE (>=),GT (>),LE (<=),LT ()
uSqr = u(x,y,0)*u(x,y,0)+u(x,y,1)*u(x,y,1)
uxy = u(x,y,0)*dble(e(i,0)) + u(x,y,1)*dble(e(i,1))
fEqei=F_eq(w(i),rhoA,uxy,uSqr)
! uxy = u(x,y,0)*dble(-e(i,0)) + u(x,y,1)*dble(-e(i,1))
uxy = -uxy ! i ，
fEqinv=F_eq(w(inv(i)),rhoA,uxy,uSqr)
f(x,y,i)=fEqei+fEqinv-fpost(x,y,inv(i))
if(f(x,y,i).lt.0.d0) f(x,y,i)=0.d0 !LT (좋은 웹페이지 즐겨찾기
