Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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