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
| !used for Loop BC treatment
subroutine cal_BC_FH(i,j, rate, f1_i,r_ff_x, r_ff_y,r_s_x,
& r_s_y, minus_i)
use ComPara, only:CP_cs_sq, fun_cal_rate, omega
use Node_data_D2Q9, only: fbar, macro, D2Q9_c, D2Q9_weight
implicit none
integer :: i,j, rate, f1_i, r_ff_x, r_ff_y, r_s_x, r_s_y,
&minus_i
real*8 :: tmp, xi, u_sf, u_w, f_star, c_squ, u_f, u_squ
u_w=0.d0
c_squ = CP_cs_sq
tmp = fun_cal_rate(i,j,rate)
! 0<q<0.5
if(tmp<0.5d0 .AND. tmp>0.d0) then
xi=omega*(2.d0*tmp-1.d0)/(1.d0-2.d0*omega)
! here is u_sf*c_i instead of u_sf
! i=1
u_sf=(D2Q9_c(f1_i,1)*macro(r_ff_x,r_ff_y,2)
&+D2Q9_c(f1_i,2)*macro(r_ff_x,r_ff_y,3))
! 0.5<=q<=1
else if (tmp>=0.5d0 .AND. tmp<1.d0) then
xi=2.d0*omega*(2.d0*tmp-1.d0)/(2.d0+omega)
u_sf=(1.d0-1.5d0/tmp)*(D2Q9_c(f1_i,1)*macro(i,j,2)
&+D2Q9_c(f1_i,2)*macro(i,j,3))+1.5d0/tmp*u_w
else
write(*,*) 'check pNode rate(k)',i,j,rate,tmp
end if
u_f=D2Q9_c(f1_i,1)*macro(i,j,2)+D2Q9_c(f1_i,2)*macro(i,j,3)
u_squ=(macro(i,j,2) ** 2.d0+macro(i,j,3)**2.d0)
f_star=D2Q9_weight(f1_i)* macro(i,j,1) * (1.d0 +
&u_sf / c_squ+ u_f ** 2.d0 / (2.d0 * c_squ ** 2.d0)
&- u_squ / (2.d0 * c_squ))
fbar(i,j,minus_i)=(1.d0-xi)*fbar(i,j,f1_i)+xi*f_star
!pNode(r_s_x,r_s_y)%node(minus_i)=xi*f_star+
!&(1.d0-xi)*pNode(i,j)%node(f1_i)
end subroutine cal_BC_FH |
