File: RELATIVE:/../../../mfix.git/model/radial_vel_correction.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: RADIAL_VEL_CORRECTION                                   !
4     !  Author: Akhilesh Bakshi                            Date: 28-JUL-14  !
5     !                                                                      !
6     !  Purpose: This module controls the iterations for solving equations  !
7     !                                                                      !
8     !  Reference: A. Bakshi, C. Altantzis, A.F. Ghoniem, "Towards accurate !
9     !  three-dimensional simulation of dense multi-phase flows using       !
10     !  cylindrical coordinates," Powder Technology, Volume 264, 09/2014,   !
11     !  Pages 242-255.                                                      !
12     !                                                                      !
13     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
14           SUBROUTINE RADIAL_VEL_CORRECTION
15     
16           use physprop
17           use fldvar
18           use constant
19           use mpi_utility
20           use functions
21     
22           IMPLICIT NONE
23     
24     ! phase index
25           INTEGER :: M
26     
27           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: Wg_Temp, Wg_Temp_GL
28           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: Ws_Temp, Ws_Temp_GL
29     
30           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: SUMX_G, SUMZ_G
31           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: SUMX_S, SUMZ_S
32     
33           INTEGER :: I1, I2, J1, K1, K2
34           INTEGER :: abIJK, abIJK1, abIJK2
35     
36           DOUBLE PRECISION :: Angle_Temp, WgAvg, WsAvg
37     
38     !-----------------------------------------------
39     ! External functions
40     !-----------------------------------------------
41           DOUBLE PRECISION, EXTERNAL :: VAVG_U_G, VAVG_V_G, VAVG_W_G, &
42                                         VAVG_U_S, VAVG_V_S, VAVG_W_S
43     
44           IF(NO_K .OR. KMAX==1) RETURN
45     
46           I1=1
47           I2=2
48     
49     ! Fukagata Scheme
50     
51     ! Implement for gas
52           ALLOCATE (Wg_Temp(ijkstart3:ijkend3))
53     
54           DO J1 = JSTART3, JEND3
55              DO K1 = KSTART3, KEND3
56              IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
57     
58                 abIJK1=FUNIJK(I2,J1,K1)
59                 Wg_Temp (abIJK1) = W_g(abIJK1)
60              END DO ! K loop Ends
61           END DO ! J Loop Ends
62     
63           ALLOCATE (Wg_Temp_GL(ijkmax3))
64           CALL GATHER (Wg_Temp, Wg_Temp_GL, root)
65           CALL BCAST (Wg_Temp_GL, root)
66     
67           ALLOCATE (SUMX_G(JMIN3:JMAX3))
68           ALLOCATE (SUMZ_G(JMIN3:JMAX3))
69     
70           DO J1=1, JMAX3
71              SUMX_G(J1)=0
72              SUMZ_G(J1)=0
73              DO K1=1, KMAX
74                 Angle_Temp = (K1-1)*2*(Pi/KMAX)
75                 IF (K1 .LE. 0.5*KMAX) THEN
76                    K2=(K1-1)+(0.5*KMAX2)
77                 ELSE
78                    K2=(K1+1)-(0.5*KMAX2)
79                 ENDIF
80     
81                 abIJK1=FUNIJK_GL (I2,J1,K1)
82                 abIJK2=FUNIJK_GL (I2,J1,K2)
83                 WgAvg=0.5*(Wg_Temp_GL(abIJK1)-Wg_Temp_GL(abIJK2))
84                 SUMX_G(J1)=SUMX_G(J1)-((2*WgAvg*sin(Angle_Temp))/KMAX)
85                 SUMZ_G(J1)=SUMZ_G(J1)+((2*WgAvg*cos(Angle_Temp))/KMAX)
86              ENDDO
87           ENDDO
88     
89           DO J1= JSTART3, JEND3
90              DO K1= KSTART3, KEND3
91                 !abIJK in local
92                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
93                 abIJK = FUNIJK (I1,J1,K1)
94                 Angle_Temp = (-Pi/KMAX)+((K1-1)*2*(Pi/KMAX))
95                 U_g(abIJK)=(SUMX_G(J1)*cos(Angle_Temp)) + &
96                     (SUMZ_G(J1)*sin(Angle_Temp))
97              ENDDO
98           ENDDO
99     
100           DEALLOCATE (Wg_Temp)
101           DEALLOCATE (Wg_Temp_GL)
102           DEALLOCATE (SUMX_G)
103           DEALLOCATE (SUMZ_G)
104     
105     
106     
107     ! Implement for Solid
108           ALLOCATE (Ws_Temp(ijkstart3:ijkend3))
109           ALLOCATE (Ws_Temp_GL(ijkmax3))
110           ALLOCATE (SUMX_S(JMIN3:JMAX3))
111           ALLOCATE (SUMZ_S(JMIN3:JMAX3))
112     
113           DO M=1, MMAX
114     
115              DO J1 = JSTART3, JEND3
116                 DO K1 = KSTART3, KEND3
117                    IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
118                    abIJK1=FUNIJK(I2,J1,K1)
119                    Ws_Temp (abIJK1) = W_s(abIJK1,M)
120                 ENDDO ! K loop Ends
121              ENDDO ! J Loop Ends
122     
123              CALL GATHER (Ws_Temp, Ws_Temp_GL, root)
124              CALL BCAST (Ws_Temp_GL, root)
125     
126              DO J1=1, JMAX3
127                 SUMX_S(J1)=0
128                 SUMZ_S(J1)=0
129                 DO K1=1, KMAX
130                    Angle_Temp = (K1-1)*2*(Pi/KMAX)
131                    IF (K1 .LE. 0.5*KMAX) THEN
132                       K2=(K1-1)+(0.5*KMAX2)
133                    ELSE
134                       K2=(K1+1)-(0.5*KMAX2)
135                    END IF
136                    abIJK1=FUNIJK_GL (I2,J1,K1)
137                    abIJK2=FUNIJK_GL (I2,J1,K2)
138                    WsAvg=0.5*(Ws_Temp_GL(abIJK1)-Ws_Temp_GL(abIJK2))
139                    SUMX_S(J1)=SUMX_S(J1)-((2*WsAvg*sin(Angle_Temp))/KMAX)
140                    SUMZ_S(J1)=SUMZ_S(J1)+((2*WsAvg*cos(Angle_Temp))/KMAX)
141                 ENDDO
142              ENDDO
143     
144              DO J1= JSTART3, JEND3
145                 DO K1= KSTART3, KEND3
146                    IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
147                    abIJK = FUNIJK (I1,J1,K1)
148                    Angle_Temp = (-Pi/KMAX)+((K1-1)*2*(Pi/KMAX))
149                    U_s(abIJK,M)=(SUMX_S(J1)*cos(Angle_Temp)) + &
150                       (SUMZ_S(J1)*sin(Angle_Temp))
151                 ENDDO
152              ENDDO
153           ENDDO ! M Loop Ends
154     
155           DEALLOCATE (Ws_Temp)
156           DEALLOCATE (Ws_Temp_GL)
157           DEALLOCATE (SUMX_S)
158           DEALLOCATE (SUMZ_S)
159     
160           RETURN
161           END SUBROUTINE RADIAL_VEL_CORRECTION
162