-
Dave Whipp authoredDave Whipp authored
define_bc_riedel.f90 4.69 KiB
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! DEFINE_BC_SPHERE Apr. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine define_bc_riedel (nnode,kfix,kfixt,x,y,z,u,v,w,temp,vo)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine assigns the boundary condition for the Stokes sphere experiment
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use constants
use definitions
implicit none
integer nnode
integer kfix(nnode*3)
integer kfixt(nnode)
double precision x(nnode),y(nnode),z(nnode)
double precision u(nnode),v(nnode),w(nnode)
double precision temp(nnode)
type (void) vo
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer i
double precision eps,delta_x,x_0,v_0,M,N
double precision, external :: ygrec
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
delta_x=0.02
x_0=.5d0
v_0=0.01d0
M=0.3
N=0.3
eps=1.d-10
do i=1,nnode
if (x(i).lt.eps) then
kfix((i-1)*3+1)=1 ; u(i)=0.d0
kfix((i-1)*3+2)=1 ; v(i)=0.d0
kfix((i-1)*3+3)=1 ; w(i)=0.d0
endif
if (x(i).gt.1.d0-eps) then
kfix((i-1)*3+1)=1 ; u(i)=0.d0 kfix((i-1)*3+2)=1 ; v(i)=0.d0
kfix((i-1)*3+3)=1 ; w(i)=0.d0
endif
if (y(i).lt.eps) then
kfix((i-1)*3+1)=1 ; u(i)=0.d0
kfix((i-1)*3+2)=1 ; v(i)=0.d0
kfix((i-1)*3+3)=1 ; w(i)=0.d0
endif
if (y(i).gt.1.d0-eps) then
kfix((i-1)*3+1)=1 ; u(i)=0.d0
kfix((i-1)*3+2)=1 ; v(i)=0.d0
kfix((i-1)*3+3)=1 ; w(i)=0.d0
endif
if (z(i).lt.eps) then
if (x(i).le.(x_0-0.5d0*M)) then
kfix((i-1)*3+2)=1 ; v(i)=-v_0/(x_0-0.5d0*M)*x(i)*ygrec(y(i),N,0.5d0)
elseif (x(i).le.(x_0-delta_x/2.d0)) then
kfix((i-1)*3+2)=1 ; v(i)=-v_0*ygrec(y(i),N,0.5d0)
elseif (x(i).le.(x_0+delta_x/2.d0)) then
kfix((i-1)*3+2)=1 ; v(i)=v_0*2.d0/pi*atan((x(i)-x_0)/delta_x)*ygrec(y(i),N,0.5d0)
elseif (x(i).le.(x_0+0.5d0*M)) then
kfix((i-1)*3+2)=1 ; v(i)=v_0*ygrec(y(i),N,0.5d0)
else
kfix((i-1)*3+2)=1 ; v(i)=v_0/(x_0-0.5d0*M)*(1.d0-x(i))*ygrec(y(i),N,0.5d0)
end if
kfix((i-1)*3+1)=1 ; u(i)=0.d0
kfix((i-1)*3+3)=1 ; w(i)=0.d0
kfixt(i)=1
temp(i)=1.d0
endif
if (z(i).gt.1.d0-eps) then
! kfix((i-1)*3+1)=1 ; u(i)=0.d0
! kfix((i-1)*3+2)=1 ; v(i)=0.d0
! kfix((i-1)*3+3)=1 ; w(i)=0.d0
kfixt(i)=1
temp(i)=0.d0
endif
if (.not.vo%influid(i)) then
kfixt(i)=1
temp(i)=0.d0
endif
end do
return
end
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
function ygrec(x,N,x0)
implicit none
double precision ygrec,x,N,x0
if (x.le.(x0-0.5d0*N)) then
ygrec=1.d0/(x0-0.5d0*N)*x
elseif (x.le.(x0+0.5d0*N)) then
ygrec=1.d0
else
ygrec=(1.d0-x)/(x0-0.5d0*N)
end if
return
end function