File: /nfs/home/0/users/jenkins/mfix.git/model/des/des_granular_temperature.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 SUBROUTINE DES_GRANULAR_TEMPERATURE
16
17
18
19
20 USE discretelement
21 USE param
22 USE param1
23 USE run
24 USE geometry
25 USE indices
26 USE compar
27 USE sendrecv
28 USE physprop
29 USE fun_avg
30 USE functions
31 USE mpi_utility
32 IMPLICIT NONE
33
34
35
36
37 INTEGER :: I, J, K, IJK
38
39 INTEGER :: M, LL
40
41 INTEGER :: NP_PHASE(DIMENSION_3, DES_MMAX)
42
43 DOUBLE PRECISION :: TEMP(DIMENSION_3, DES_MMAX)
44
45 INTEGER :: PC
46
47 DOUBLE PRECISION :: SQR_VEL, SQR_ROT_VEL
48
49
50
51
52
53
54
55
56
57
58 (:,:) = ZERO
59 NP_PHASE(:,:) = ZERO
60 PC = 0
61 DO LL = 1, MAX_PIP
62
63 IF(.NOT.PEA(LL,1)) CYCLE
64 PC = PC + 1
65
66 IF(PEA(LL,4)) CYCLE
67
68 I = PIJK(LL,1)
69 J = PIJK(LL,2)
70 K = PIJK(LL,3)
71 IJK = PIJK(LL,4)
72
73 M = PIJK(LL,5)
74 NP_PHASE(IJK,M) = NP_PHASE(IJK,M) + 1
75
76 TEMP(IJK,M) = TEMP(IJK,M) + &
77 (DES_VEL_NEW(1,LL)-DES_U_s(IJK,M))**2
78 TEMP(IJK,M) = TEMP(IJK,M) + &
79 (DES_VEL_NEW(2,LL)-DES_V_s(IJK,M))**2
80 IF(DO_K) THEN
81 TEMP(IJK,M) = TEMP(IJK,M) + &
82 (DES_VEL_NEW(3,LL)-DES_W_s(IJK,M))**2
83 ENDIF
84
85 IF(PC .EQ. PIP) EXIT
86 ENDDO
87
88
89 DO IJK = IJKSTART3, IJKEND3
90 IF(FLUID_AT(IJK)) THEN
91
92 DO M = 1,DES_MMAX
93 IF (NP_PHASE(IJK,M) > 0 ) THEN
94 DES_THETA(IJK,M) = TEMP(IJK,M)/&
95 DBLE(DIMN*NP_PHASE(IJK,M))
96 ELSE
97 DES_THETA(IJK,M) = ZERO
98 ENDIF
99 ENDDO
100 ENDIF
101 ENDDO
102
103
104
105
106
107
108
109 = ZERO
110 DES_ROTE = ZERO
111 DES_PE = ZERO
112 DES_VEL_AVG(:) = ZERO
113
114
115
116 = 0
117 DO LL = 1, MAX_PIP
118 IF(PEA(LL,1)) PC = PC + 1
119
120 IF(.NOT.PEA(LL,1) .OR. PEA(LL,4)) CYCLE
121
122 SQR_VEL = ZERO
123 SQR_ROT_VEL = ZERO
124 DO I = 1, DIMN
125 SQR_VEL = SQR_VEL + DES_VEL_NEW(I,LL)**2
126 ENDDO
127
128 DO I = 1, merge(1,3,NO_K)
129 SQR_ROT_VEL = SQR_ROT_VEL + OMEGA_NEW(I,LL)**2
130 ENDDO
131
132
133 DES_KE = DES_KE + PMASS(LL)/2.d0 * SQR_VEL
134
135 = DES_ROTE + &
136 (0.4D0*PMASS(LL)*DES_RADIUS(LL)**2)/2.d0 * SQR_ROT_VEL
137 DES_PE = DES_PE + PMASS(LL)*DBLE(ABS(GRAV(2)))*&
138 DES_POS_NEW(2,LL)
139 DES_VEL_AVG(:) = DES_VEL_AVG(:) + DES_VEL_NEW(:,LL)
140
141 IF(PC .EQ. PIP) EXIT
142 ENDDO
143
144
145
146 CALL GLOBAL_ALL_SUM(PIP-iGHOST_CNT,TOT_PAR)
147 CALL GLOBAL_ALL_SUM(DES_VEL_AVG(1:DIMN))
148
149
150 IF(TOT_PAR > 0) DES_VEL_AVG(:) = DES_VEL_AVG(:)/DBLE(TOT_PAR)
151
152
153
154
155
156 (:) = ZERO
157 GLOBAL_GRAN_TEMP(:) = ZERO
158 PC = 0
159 DO LL = 1, MAX_PIP
160 IF(PEA(LL,1)) PC = PC + 1
161
162 IF(.NOT.PEA(LL,1) .OR. PEA(LL,4)) CYCLE
163
164 GLOBAL_GRAN_ENERGY(:) = GLOBAL_GRAN_ENERGY(:) + &
165 0.5d0*PMASS(LL)*(DES_VEL_NEW(:,LL)-DES_VEL_AVG(:))**2
166 GLOBAL_GRAN_TEMP(:) = GLOBAL_GRAN_TEMP(:) + &
167 (DES_VEL_NEW(:,LL)-DES_VEL_AVG(:))**2
168
169 IF(PC .EQ. PIP) EXIT
170 ENDDO
171
172 CALL GLOBAL_ALL_SUM(GLOBAL_GRAN_TEMP)
173 CALL GLOBAL_ALL_SUM(GLOBAL_GRAN_ENERGY)
174 CALL GLOBAL_ALL_SUM(DES_KE)
175 CALL GLOBAL_ALL_SUM(DES_PE)
176 CALL GLOBAL_ALL_SUM(DES_ROTE)
177
178 IF(TOT_PAR > 0) GLOBAL_GRAN_ENERGY(:) = &
179 GLOBAL_GRAN_ENERGY(:)/DBLE(TOT_PAR)
180 IF(TOT_PAR > 0) GLOBAL_GRAN_TEMP(:) = &
181 GLOBAL_GRAN_TEMP(:)/DBLE(TOT_PAR)
182
183
184
185 RETURN
186
187 END SUBROUTINE DES_GRANULAR_TEMPERATURE
188
189
190
191
192
193
194
195
196
197
198
199 SUBROUTINE CALC_DES_BEDHEIGHT
200
201
202
203
204 USE indices
205 USE geometry
206 USE compar
207 USE discretelement
208 USE des_bc
209 USE physprop
210 USE fldvar
211 USE functions
212 IMPLICIT NONE
213
214
215
216
217 INTEGER :: L
218
219 INTEGER :: PC
220
221 INTEGER :: M
222
223 INTEGER :: J, IJK
224
225 DOUBLE PRECISION :: hcell
226
227 DOUBLE PRECISION :: hpart
228
229 DOUBLE PRECISION :: EP_SM
230
231 DOUBLE PRECISION :: tmp_num(DES_MMAX), tmp_den(DES_MMAX)
232
233
234
235
236 (:) = ZERO
237 tmp_den(:) = ZERO
238 DO IJK = IJKSTART3, IJKEND3
239 J = J_OF(IJK)
240 DO M = 1, DES_MMAX
241 IF(DES_ROP_S(IJK,M) > ZERO) THEN
242 hcell = 0.5d0*(YN(J)+YN(J-1))
243 EP_SM = DES_ROP_S(IJK,M)/DES_RO_S(M)
244 tmp_num(M) = tmp_num(M) + EP_SM*hcell*VOL(IJK)
245 tmp_den(M) = tmp_den(M) + EP_SM*VOL(IJK)
246 ENDIF
247 ENDDO
248 ENDDO
249
250 IF (PIP >0) bed_height(:) = tmp_num(:)/tmp_den(:)
251
252
253
254 IF(.FALSE.) THEN
255 tmp_num(:) = ZERO
256 tmp_den(:) = ZERO
257
258 PC = 0
259 DO L = 1, MAX_PIP
260 IF(PEA(L,1)) PC = PC + 1
261
262 IF(.NOT.PEA(L,1) .OR. PEA(L,4)) CYCLE
263 M = PIJK(L,5)
264 hpart = DES_POS_NEW(2,L)
265 tmp_num(M) = tmp_num(M) + hpart
266 tmp_den(M) = tmp_den(M) + 1
267 IF(PC .EQ. PIP) EXIT
268 ENDDO
269
270 IF (PIP >0) bed_height(:) = tmp_num(:)/tmp_den(:)
271 ENDIF
272
273 RETURN
274
275 END SUBROUTINE CALC_DES_BEDHEIGHT
276
277