cavity.f90

Fortran语言

module const
implicit none
real length_x,length_y,dx,dy,re,u0,tt,gravity
integer mx,my
end module const

module variables
implicit none
real dt
type type_var
	real u,v,p,xi,eta
end type
type (type_var),allocatable :: fluid(:,:)
end module variables

!*******************Main program **************************
program main
use const
use variables
implicit none
real tempt
integer nn
tempt=0.0
nn=0
call input

call initial

call output(0)

do while(tt.gt.tempt)
	nn=nn+1

!	call detat
	dt = 0.001
	tempt=tempt+dt
	write(*,*)"Step:",nn,"  dt=",dt,"   time=",tempt

	call pressure
 	

	call velocity

	call impressibility
	
	if(mod(nn,100).eq.0) call output(nn)


end do
write(*,*) "***************** OK ********************"
end program main
!*********************End main program ********************


!******************* Subroutine input**********************
subroutine input
use const
implicit none
open(2,file="input.dat",status="old")

read(2,*) length_x

read(2,*) length_y

read(2,*) u0

read(2,*) re

read(2,*) mx,my

dx=length_x/mx
dy=length_y/my

read(2,*) tt

close(2)
end subroutine input
!****************End subroutine input *********************


!****************** Subroutine initial ********************
subroutine initial
use const
use variables
implicit none
integer i,j,allocate_error

allocate(fluid(mx+3,my+3),stat=allocate_error) 

do j=1,my+3
do i=1,mx+3
	fluid(i,j)%u=0.0
	fluid(i,j)%v=0.0
	fluid(i,j)%p=0.0
end do
end do

do i=3,mx+1
	fluid(i,my+1)%u=0.5*u0
	fluid(i,my+2)%u=1.5*u0
end do

gravity = 0.0

call bound_velocity

call bound_press

end subroutine initial
!****************** End subroutine initial **************** 


!****************** Subroutine Bound_velocity *************
subroutine bound_velocity
use const
use variables
implicit none
integer i,j

do i=1,mx+3
  !******lower
	fluid(i,2)%v=0.0
	fluid(i,1)%v= fluid(i,3)%v
	fluid(i,1)%u=-fluid(i,2)%u
  !******upper
    fluid(i,my+2)%v=0.0
	fluid(i,my+3)%v= fluid(i,my+1)%v
	fluid(i,my+2)%u=-fluid(i,my+1)%u+2*u0
end do

do j=1,my+3
  !******left
	fluid(2,j)%u=0.0
	fluid(1,j)%u= fluid(3,j)%u
	fluid(1,j)%v=-fluid(2,j)%v
  !******right
	fluid(mx+2,j)%u=0.0
	fluid(mx+3,j)%u= fluid(mx+1,j)%u
	fluid(mx+2,j)%v=-fluid(mx+1,j)%v
end do


end subroutine bound_velocity
!****************** End subroutine Bound_velocity *********

!****************** Subroutine detat **********************
subroutine detat
use const
use variables
implicit none

real dt1,dt2,cfl1,cfl2
real temp1,temp2
integer i,j,ii,jj
temp1=0.0
cfl1=0.2
cfl2=0.5

do j=1,my+2
do i=1,mx+2
	temp2=abs(fluid(i,j)%u)/dx+abs(fluid(i,j)%v)/dy
	if(temp2.gt.temp1) 	then
		temp1=temp2
		ii=i
		jj=j
	end if	
end do
end do

dt1=cfl1/temp1
dt2=cfl2*re*0.5*(dx**2)*(dy**2)/(dx**2+dy**2)

!write(*,*) "    min deta ",fluid(ii,jj)%u,fluid(ii,jj)%v,ii,jj
dt=min(dt1,dt2)
end subroutine detat
!********************* End subroutine detat ***************


!******************** Subroutine xi_eta ***************
subroutine xi_eta
use const
use variables
implicit none
integer i,j
real temp,temp1

do j=2,my+1
do i=2,mx+2
	temp=((fluid(i-1,j)%u-2*fluid(i,j)%u+fluid(i+1,j)%u)/dx**2+ &
		  (fluid(i,j-1)%u-2*fluid(i,j)%u+fluid(i,j+1)%u)/dy**2)/re
	temp1=((fluid(i,j)%u+fluid(i+1,j)%u)**2-(fluid(i-1,j)%u+fluid(i,j)%u)**2)/dx
	temp1=temp1+((fluid(i,j)%u+fluid(i,j+1)%u)*(fluid(i-1,j+1)%v+fluid(i,j+1)%v)-	    &
			     (fluid(i,j)%u+fluid(i,j-1)%u)*(fluid(i-1,j  )%v+fluid(i,j  )%v))/dy
	fluid(i,j)%xi=temp-0.25*temp1	
end do
end do

do j=2,my+2
do i=2,mx+1
	temp=((fluid(i-1,j)%v-2*fluid(i,j)%v+fluid(i+1,j)%v)/dx**2+ &
		  (fluid(i,j-1)%v-2*fluid(i,j)%v+fluid(i,j+1)%v)/dy**2)/re
	temp1=((fluid(i,j)%v+fluid(i,j+1)%v)**2-(fluid(i,j-1)%v+fluid(i,j)%v)**2)/dy
	temp1=temp1+((fluid(i+1,j)%u+fluid(i+1,j-1)%u)*(fluid(i,j)%v+fluid(i+1,j)%v)-    &
				 (fluid(i  ,j)%u+fluid(i  ,j-1)%u)*(fluid(i,j)%v+fluid(i-1,j)%v))/dx
	fluid(i,j)%eta=temp-0.25*temp1-gravity
end do
end do

end subroutine xi_eta
!******************** End subroutine xi_eta ***********


!******************** Subroutine pressure *****************
subroutine pressure
use const
use variables
implicit none
integer i,j,iteration
real pmin,wopt,temp,error,max,x,y

real,allocatable::rhs(:,:)

allocate(rhs(2:mx+1,2:my+1))

call xi_eta

do j=2,my+1
do i=2,mx+1
   rhs(i,j)=(fluid(i+1,j)%xi-fluid(i,j)%xi) /dx+(fluid(i,j+1)%eta-fluid(i,j)%eta)/dy &
		   +(fluid(i+1,j)%u-fluid(i,j)%u)/dt/dx+(fluid(i,j+1)%v-fluid(i,j)%v)/dt/dy	
enddo
enddo

open(22,file='rhs.dat')
	write(22,*) ' Title="The cavity flow" '
	write(22,*) ' Variables="x" "y" "press" '
	write(22,*) ' Zone I=',mx,' , J=',my,'  , F=Point'
	do j=2,my+1
	do i=2,mx+1
		x=(i-1.5)*dx
		y=(j-1.5)*dy
		write(22,*) x,y,rhs(i,j)
	enddo
	enddo
close(22)

pmin=1.0e-7
wopt=1.93
error=1000
iteration=0
max=25000
do while(error.gt.pmin)
	iteration=iteration+1
	error=0
	call bound_press
	do j=2,my+1
	do i=2,mx+1
		temp=(fluid(i+1,j)%p+fluid(i-1,j)%p)/dx**2+(fluid(i,j+1)%p+fluid(i,j-1)%p)/dy**2
		temp=(temp-rhs(i,j))/(2/dx**2+2/dy**2)
		temp=wopt*(temp-fluid(i,j)%p)
		if(abs(temp).gt.error) error=abs(temp)
		fluid(i,j)%p=temp+fluid(i,j)%p
	end do
	end do
	if(mod(iteration,1000).eq.0) then
		write(*,*) "    Iteration:",iteration,"    error:",error
	end if
	if(iteration.gt.max) stop "No convergence"
end do

write(*,*) "iteration:  ",iteration,"  error:",error

deallocate(rhs)
end subroutine pressure
!******************** End subroutine pressure *************


!******************** Subroutine bound_press  *************
subroutine bound_press
use const
use variables
implicit none
integer i,j
real temp

do j=2,my+1
	temp=((fluid(1,j  )%u-2*fluid(2,j)%u+fluid(3,j  )%u)/dx**2+ &
		  (fluid(2,j-1)%u-2*fluid(2,j)%u+fluid(2,j+1)%u)/dy**2)/re
	fluid(1,j)%p=fluid(2,j)%p-dx*temp
	
	temp=((fluid(mx+1,j  )%u-2*fluid(mx+2,j)%u+fluid(mx+3,j  )%u)/dx**2+ &
		  (fluid(mx+2,j-1)%u-2*fluid(mx+2,j)%u+fluid(mx+2,j+1)%u)/dy**2)/re
	fluid(mx+2,j)%p=fluid(mx+1,j)%p+dx*temp	
end do

do i=2,mx+1
	temp=((fluid(i-1,2)%v-2*fluid(i,2)%v+fluid(i+1,2)%v)/dx**2+ &
		  (fluid(i  ,1)%v-2*fluid(i,2)%v+fluid(i  ,3)%v)/dy**2)/re-gravity
	fluid(i,1)%p=fluid(i,2)%p-dy*temp
	
	temp=((fluid(i-1,my+2)%v-2*fluid(i,my+2)%v+fluid(i+1,my+2)%v)/dx**2+ &
		  (fluid(i  ,my+1)%v-2*fluid(i,my+2)%v+fluid(i  ,my+3)%v)/dy**2)/re-gravity
	fluid(i,my+2)%p=fluid(i,my+1)%p+dy*temp	
end do

end subroutine bound_press
!******************** End subroutine bound_press **********

!******************** Subroutine Velocity *****************
subroutine velocity
use const
use variables
implicit none
integer i,j
real temp

do j=2,my+1
do i=3,mx+1
	temp=dt*(-(fluid(i,j)%p-fluid(i-1,j)%p)/dx+fluid(i,j)%xi)
	fluid(i,j)%u=fluid(i,j)%u+temp
end do
end do

do j=3,my+1
do i=2,mx+1
	temp=dt*(-(fluid(i,j)%p-fluid(i,j-1)%p)/dy+fluid(i,j)%eta)
	fluid(i,j)%v=fluid(i,j)%v+temp
end do
end do

call bound_velocity

end subroutine velocity
!******************** End subroutine velociy **************

!******************** Subroutine output *******************
subroutine output(nn)
use const
use variables
implicit none
intent(in) nn
integer i,j,nn
real x,y
character*5 name

write(name,'(i5.5)') nn

open(2,file="dat"//name//'.dat')
write(2,*) ' Title="The cavity flow" '
write(2,*) ' Variables="x" "y" "press" "u" "v" "xci" "eta" "px" "py" '
write(2,*) ' Zone I=',mx,' , J=',my,'  , F=Point'
do j=2,my+1
do i=2,mx+1
	x=(i-1.5)*dx
	y=(j-1.5)*dy
	write(2,*) x,y,fluid(i,j)%p,(fluid(i,j)%u+fluid(i+1,j)%u)/2 &
							   ,(fluid(i,j)%v+fluid(i,j+1)%v)/2	&
							   ,fluid(i,j)%xi,fluid(i,j)%eta	&
							   ,(fluid(i,j)%p-fluid(i-1,j)%p)/dx &
							   ,(fluid(i,j)%p-fluid(i,j-1)%p)/dy
end do 
end do
close(2)

end subroutine output
!******************** End subroutine output ***************


!******************** Subroutine impressibility ***********
subroutine impressibility
use const
use variables
implicit none
integer i,j,ii,jj
real error,temp

error=0
do j=2,my+1
do i=2,mx+1
	temp=(fluid(i+1,j)%u-fluid(i,j)%u)/dx+(fluid(i,j+1)%v-fluid(i,j)%v)/dy
	if(abs(temp).gt.abs(error)) error=temp
	ii=i
	jj=j
end do
end do
write(*,*) "    impressibility:",error,ii,jj
end subroutine impressibility
!******************** End subroutine impressibility *******