File: /nfs/home/0/users/jenkins/mfix.git/model/des/set_bc_dem_mi.f
1
2
3
4
5
6
7
8
9 SUBROUTINE SET_BC_DEM_MI
10
11
12
13
14 USE bc
15 USE compar
16 USE constant
17 USE des_allocate
18 USE des_bc
19 USE discretelement
20 USE error_manager
21 USE funits
22 USE geometry
23 USE indices
24 USE param
25 USE param1
26 USE physprop
27 USE run
28 USE toleranc
29
30 IMPLICIT NONE
31
32 INTEGER :: BCV
33
34
35
36
37 INTEGER BCV_I
38 INTEGER M, MM
39 INTEGER HOLD, I
40 INTEGER RANGE_TOP, RANGE_BOT
41 INTEGER PHASE_CNT
42 INTEGER PHASE_LIST(DIM_M)
43
44
45 DOUBLE PRECISION NPMpSEC(DES_MMAX)
46 DOUBLE PRECISION NPpSEC
47 DOUBLE PRECISION NPpDT
48 DOUBLE PRECISION SCALED_VAL
49 DOUBLE PRECISION MAX_DIA
50
51 DOUBLE PRECISION :: EPs_ERR
52 DOUBLE PRECISION :: VOL_FLOW
53
54 LOGICAL, parameter :: setDBG = .FALSE.
55 LOGICAL :: dFlag
56
57 LOGICAL :: FATAL
58
59
60 DOUBLE PRECISION VEL_TMP(DIM_M)
61 DOUBLE PRECISION EPs_TMP(DIM_M)
62
63
64
65 DOUBLE PRECISION MIN_VEL, MAX_VEL
66
67 DOUBLE PRECISION MINIPV, MAXIPV
68
69 INTEGER :: OCCUPANTS
70
71
72
73
74 CALL INIT_ERR_MSG("SET_BC_DEM_MI")
75
76 dFlag = (DMP_LOG .AND. setDBG)
77 if(dFlag) write(*,"(2/,2x,'DEM inlet count: ',I4)") DEM_BCMI
78
79
80
81 IF(RUN_TYPE /= 'RESTART_1') CALL ALLOCATE_DEM_MI
82
83
84 DO BCV_I = 1, DEM_BCMI
85
86
87 = DEM_BCMI_MAP(BCV_I)
88
89 if(dFlag) write(*,"(2/,'Setting DEM_MI:',I3)") BCV_I
90
91
92
93 = 0
94
95 (:) = -1
96
97 = ZERO
98
99
100 DO M=1, SMAX + DES_MMAX
101 IF(SOLIDS_MODEL(M) /= 'DEM') CYCLE
102 IF(BC_ROP_s(BCV,M) == UNDEFINED) CYCLE
103 IF(COMPARE(BC_ROP_s(BCV,M),ZERO)) CYCLE
104 PHASE_CNT = PHASE_CNT + 1
105 PHASE_LIST(PHASE_CNT) = M-SMAX
106 MAX_DIA = MAX(MAX_DIA,DES_D_P0(M-SMAX))
107 ENDDO
108
109 (BCV_I)%POLYDISPERSE = (PHASE_CNT > 1)
110
111
112 CALL LAYOUT_MI_DEM(BCV, BCV_I, MAX_DIA)
113
114
115 = ZERO
116 NPMpSEC(:) = ZERO
117
118
119 DO MM = 1, PHASE_CNT
120 M = PHASE_LIST(MM)
121
122
123 SELECT CASE(BC_PLANE(BCV))
124 CASE('N','S'); VEL_TMP(M) = abs(BC_V_s(BCV,M+SMAX))
125 CASE('E','W'); VEL_TMP(M) = abs(BC_U_s(BCV,M+SMAX))
126 CASE('T','B'); VEL_TMP(M) = abs(BC_W_s(BCV,M+SMAX))
127 END SELECT
128
129
130 = MAX(ABS(VEL_TMP(M)), MAX_VEL)
131
132
133
134 = VEL_TMP(M) * BC_AREA(BCV) * BC_EP_S(BCV,M+SMAX)
135
136
137 (M) = VOL_FLOW / (PI/6.d0*DES_D_P0(M)**3)
138
139 if(dFlag) write(*,1100) M, VEL_TMP(M), NPMpSEC(M)
140 ENDDO
141
142
143 = sum(NPMpSEC)
144
145 1100 FORMAT(/2x,'Conversion Info: Phase ',I2,/4x,'Velocity: ',g12.5,/ &
146 4X,'NPMpSEC = ',F11.1)
147
148 if(dFlag) write(*,"(/2x,'Max Velocity:',3x,g12.5)") MAX_VEL
149 if(dFlag) write(*,"( 2x,'NPpSEC:',3x,F11.1)") NPpSEC
150
151
152 = NPpSEC * DTSOLID
153
154
155 IF(NPpDT .LT. 1)THEN
156 PI_COUNT(BCV_I) = 1
157 PI_FACTOR(BCV_I) = FLOOR(real(1.d0/NPpDT))
158
159 ELSE
160 PI_COUNT(BCV_I) = CEILING(real(NPpDT))
161 PI_FACTOR(BCV_I) = 1
162 ENDIF
163
164 OCCUPANTS = DEM_MI(BCV_I)%OCCUPANTS
165
166
167
168 = MAX_DIA/( DTSOLID*dble(PI_FACTOR(BCV_I)) * dble( &
169 FLOOR( real(OCCUPANTS)/real(PI_COUNT(BCV_I)))))
170
171
172
173 = MAX_DIA/( DTSOLID*dble(PI_FACTOR(BCV_I)) * dble( &
174 FLOOR(CEILING(real(OCCUPANTS)/2.0)/real(PI_COUNT(BCV_I)))))
175
176 if(dFlag) write(*,"(/2x,'MaxIPV:',3x,g12.5)") MAXIPV
177 if(dFlag) write(*,"( 2x,'MinIPV:',3x,g12.5)") MINIPV
178
179 IF(MAX_VEL < MINIPV) THEN
180 WRITE(ERR_MSG,1110) BCV, MAX_VEL, MINIPV
181 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
182 ENDIF
183
184 1110 FORMAT('Error 1110: Solids velocity for BC ',I3,' is too low ', &
185 'to satisfy DEM',/'inlet packing restrictions. Potential ', &
186 'solutions:',//,' > If the BC velocities (BC_U_s, BC_V_s, ', &
187 'BC_W_s) are defined, specify',/3x,'a larger value for the ',&
188 'velocity normal to the flow plane.',//' > If MASSFLOW or ', &
189 'VOLFLOW are defined, decrease the solids volume',/3x, &
190 'fraction to increase solids velocity.',//2x,'Max user-', &
191 'specified BC velocity: ',g12.5,/2x,'Minimum required ', &
192 'solids Velocity: ',g12.5)
193
194
195
196 = ZERO
197 EPs_TMP = ZERO
198 DO MM = 1, PHASE_CNT
199 M = PHASE_LIST(MM)
200 EPs_TMP(M) = BC_EP_s(BCV,M+SMAX) * (VEL_TMP(M) / MAX_VEL)
201 EPs_ERR = EPs_ERR + (BC_EP_s(BCV,M+SMAX) - EPs_TMP(M))
202
203
204 SELECT CASE(BC_PLANE(BCV))
205 CASE('N'); BC_V_s(BCV,M) = abs(MAX_VEL)
206 CASE('S'); BC_V_s(BCV,M) = -abs(MAX_VEL)
207 CASE('E'); BC_U_s(BCV,M) = abs(MAX_VEL)
208 CASE('W'); BC_U_s(BCV,M) = -abs(MAX_VEL)
209 CASE('T'); BC_W_s(BCV,M) = abs(MAX_VEL)
210 CASE('B'); BC_W_s(BCV,M) = -abs(MAX_VEL)
211 END SELECT
212
213 ENDDO
214
215
216
217 IF(.NOT.COMPARE(EPs_ERR,SMALL_NUMBER)) THEN
218 IF(EPs_ERR > 0.01) THEN
219 WRITE(ERR_MSG,1200) BCV, MAX_VEL, EPs_ERR
220 CALL FLUSH_ERR_MSG(FOOTER=.FALSE.)
221 FATAL = .TRUE.
222
223
224
225 ELSE
226 WRITE(ERR_MSG,1205) BCV, MAX_VEL, EPs_ERR
227 CALL FLUSH_ERR_MSG(FOOTER=.FALSE.)
228 FATAL = .FALSE.
229 ENDIF
230
231 WRITE(ERR_MSG, 1210)
232 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
233 DO MM = 1, PHASE_CNT
234 M = PHASE_LIST(MM)
235 WRITE(ERR_MSG,1211) M, BC_EP_s(BCV,M), EPs_TMP(M), &
236 (BC_EP_s(BCV,M)-EPs_TMP(M))
237 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
238 ENDDO
239 CALL FLUSH_ERR_MSG(HEADER=.FALSE., ABORT=FATAL)
240 ENDIF
241
242
243 1200 FORMAT('Error 1200: Unable to impose a uniform solids velocity ',&
244 'on BC',I3,'.',/'Setting all solids to the highest velocity ',&
245 'results in a large error',/'in the solids volume fraction. ',&
246 'Please correct the mfix.dat file.',2/5X,'Max Inlet Velocity',&
247 ':',1X,ES11.4,/5x,'Total BC_EP_s Error:',ES11.4)
248
249
250 1205 FORMAT('Warning 1201: Uniform solids velocity imposed on BC',I3, &
251 '.',2/,5X,'Uniform Inlet Velocity:',1X,ES11.4,/5X,'Total ', &
252 'BC_EP_s Error:',4X,ES11.4,/' ')
253
254 1210 FORMAT(/,5X,'|',11('-'),'|',3(14('-'),'|'),/5X,'|',3X,'Phase',3X,&
255 '|',4X,'BC_EP_s',3X,'|',2X,'Calculated',2X,'|',3X,'ABS ', &
256 'Error',2X,'|',/5X,'|',11('-'),'|',3(14('-'),'|'))
257
258 1211 FORMAT(5X,'|',4X,I2,5X,'| ',1X,ES11.4,1X,'|',1X,ES11.4,2X,'|', &
259 1X,ES11.4,1X,' |',/5X,'|',11('-'),'|',3(14('-'),'|'))
260
261
262
263 IF(DEM_MI(BCV_I)%POLYDISPERSE) THEN
264 RANGE_BOT = 1
265 DO MM=1,PHASE_CNT - 1
266 M = PHASE_LIST(MM)
267 SCALED_VAL = dble(NUMFRAC_LIMIT)*(NPMpSEC(M)/NPpSEC)
268 RANGE_TOP = FLOOR(SCALED_VAL) + (RANGE_BOT-1)
269 DEM_BC_POLY_LAYOUT(BCV_I,RANGE_BOT:RANGE_TOP) = M
270 RANGE_BOT = RANGE_TOP+1
271 ENDDO
272
273 M = PHASE_LIST(PHASE_CNT)
274 DEM_BC_POLY_LAYOUT(BCV_I,RANGE_BOT:NUMFRAC_LIMIT) = M
275
276 ELSE
277 DEM_BC_POLY_LAYOUT(BCV_I,:) = PHASE_LIST(1)
278 ENDIF
279
280
281
282
283
284 IF(RUN_TYPE /= 'RESTART_1') &
285 DEM_MI_TIME(BCV_I) = TIME + dble(PI_FACTOR(BCV_I)) * DTSOLID
286
287
288 WRITE(ERR_MSG,1000) BCV, NPpDT, int(NPpDT/DTSOLID), &
289 PI_FACTOR(BCV_I), PI_COUNT(BCV_I), DEM_MI_TIME(BCV_I)
290 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
291
292 1000 FORMAT(2/,2X,'For mass inlet BC: ', I3,/,&
293 4X,'No. particles injected per solids time step = ', ES15.8,/,&
294 4X,'No. particles injected per second = ', I10,/,&
295 4X,'PI_FACTOR = ', I10,' PI_COUNT = ', I5,/,&
296 4X,'start DES_MI_TIME = ', ES15.8,/' ')
297
298 ENDDO
299
300
301 CALL SET_DEM_BCMI_IJK
302
303
304 CALL FINL_ERR_MSG
305
306
307 RETURN
308 END SUBROUTINE SET_BC_DEM_MI
309
310
311
312
313
314
315
316
317
318
319
320
321
322 SUBROUTINE SET_DEM_BCMI_IJK
323
324 use bc, only: BC_PLANE
325 use bc, only: BC_I_w, BC_I_e, BC_J_s, BC_J_n, BC_K_b, BC_K_t
326
327 use des_bc, only: DEM_BCMI, DEM_BCMI_MAP, DEM_BCMI_IJK
328 use des_bc, only: DEM_BCMI_IJKSTART, DEM_BCMI_IJKEND
329
330 use funits, only: DMP_LOG
331
332 use mpi_utility
333 use error_manager
334 use functions
335
336 IMPLICIT NONE
337
338 INTEGER, ALLOCATABLE :: LOC_DEM_BCMI_IJK(:)
339
340 INTEGER :: BCV, BCV_I
341
342 INTEGER :: LC
343
344 INTEGER :: MAX_CELLS
345
346 INTEGER :: BND1, BND2
347
348 LOGICAL, parameter :: setDBG = .FALSE.
349 LOGICAL :: dFlag
350
351 INTEGER :: I,J,K,IJK
352 INTEGER :: I_w, I_e, J_s, J_n, K_b, K_t
353
354 CALL INIT_ERR_MSG("SET_DEM_BCMI_IJK")
355
356 dFlag = (DMP_LOG .AND. setDBG)
357
358 if(dFlag) write(*,"(2/,2x,'From: SET_DEM_BCMI_IJK')")
359
360
361
362
363
364 = 0
365 DO BCV_I=1, DEM_BCMI
366 BCV = DEM_BCMI_MAP(BCV_I)
367
368
369 if(dFlag) WRITE(*,"(/2x,'Adding cells for BCV: ',I3)") BCV
370 SELECT CASE (BC_PLANE(BCV))
371 CASE('N','S')
372 BND1 = min(BC_I_e(BCV)+1,IMAX1) - max(BC_I_w(BCV)-1,IMIN1)
373 BND2 = min(BC_K_t(BCV)+1,KMAX1) - max(BC_K_b(BCV)-1,KMIN1)
374
375 CASE('E','W')
376 BND1 = min(BC_J_n(BCV)+1,JMAX1) - max(BC_J_s(BCV)-1,JMIN1)
377 BND2 = min(BC_K_t(BCV)+1,KMAX1) - max(BC_K_b(BCV)-1,KMIN1)
378
379 CASE('T','B')
380 BND1 = min(BC_I_e(BCV)+1,IMAX1) - max(BC_I_w(BCV)-1,IMIN1)
381 BND2 = min(BC_J_n(BCV)+1,JMAX1) - max(BC_J_s(BCV)-1,JMIN1)
382 END SELECT
383
384 MAX_CELLS = MAX_CELLS + (BND1 + 1)*(BND2 + 1)
385 if(dFlag) WRITE(*,"(4x,'Plane: ',A)") BC_PLANE(BCV)
386 if(dFlag) WRITE(*,"(4x,'Cells: ', I8)") (BND1+1)*(BND2+1)
387 ENDDO
388
389 if(dFlag) write(*,"(2x,'Max Cells: ',I8)") MAX_CELLS
390
391
392
393 allocate( LOC_DEM_BCMI_IJK(MAX_CELLS) )
394
395
396
397
398 = 1
399 DO BCV_I = 1, DEM_BCMI
400
401 DEM_BCMI_IJKSTART(BCV_I) = LC
402 BCV = DEM_BCMI_MAP(BCV_I)
403
404 if(dFlag) write(*,"(/2x,'Searching for fluid cells:',I3)") BCV
405
406 I_w = max(BC_I_w(BCV)-1,IMIN1); I_e = min(BC_I_e(BCV)+1,IMAX1)
407 J_s = max(BC_J_s(BCV)-1,JMIN1); J_n = min(BC_J_n(BCV)+1,JMAX1)
408 K_b = max(BC_K_b(BCV)-1,KMIN1); K_t = min(BC_K_t(BCV)+1,KMAX1)
409
410
411
412 SELECT CASE (BC_PLANE(BCV))
413 CASE('N'); J_s = BC_J_s(BCV)+1; J_n = BC_J_n(BCV)+1
414 CASE('S'); J_s = BC_J_s(BCV)-1; J_n = BC_J_n(BCV)-1
415 CASE('E'); I_w = BC_I_w(BCV)+1; I_e = BC_I_e(BCV)+1
416 CASE('W'); I_w = BC_I_w(BCV)-1; I_e = BC_I_e(BCV)-1
417 CASE('T'); K_b = BC_K_b(BCV)+1; K_t = BC_K_t(BCV)+1
418 CASE('B'); K_b = BC_K_b(BCV)-1; K_t = BC_K_t(BCV)-1
419 END SELECT
420
421 if(dFlag) then
422 write(*,"(4x,'Search bounds: ')")
423 write(*,"(6x,'I_w/I_e:',2(2x,I6))") I_w, I_e
424 write(*,"(6x,'J_s/J_n:',2(2x,I6))") J_s, J_n
425 write(*,"(6x,'K_b/K_t:',2(2x,I6))") K_b, K_t
426 endif
427
428
429 DO K = K_b, K_t
430 DO J = J_s, J_n
431 DO I = I_w, I_e
432
433 IF(.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
434 IF(DEAD_CELL_AT(I,J,K)) CYCLE
435
436 IJK = FUNIJK(I,J,K)
437 LOC_DEM_BCMI_IJK(LC) = IJK
438 LC = LC+1
439 ENDDO
440 ENDDO
441 ENDDO
442
443 DEM_BCMI_IJKEND(BCV_I) = LC-1
444
445 IF(dFLAG) write(*,1111) BCV, BCV_I, &
446 DEM_BCMI_IJKSTART(BCV_I), DEM_BCMI_IJKEND(BCV_I)
447
448 ENDDO
449
450 1111 FORMAT(/2x,'DEM Mass Inflow:',/4x,'BC:',I4,3x,'MAP:',I4, &
451 /4x,'IJKSTART:',I6,/4x,'IJKEND: ',I6)
452
453
454
455 IF(LC > 1) THEN
456 allocate( DEM_BCMI_IJK(LC-1) )
457 DEM_BCMI_IJK(1:LC-1) = LOC_DEM_BCMI_IJK(1:LC-1)
458 ELSE
459 allocate( DEM_BCMI_IJK(1) )
460 DEM_BCMI_IJK(1) = LOC_DEM_BCMI_IJK(1)
461 ENDIF
462
463 deallocate(LOC_DEM_BCMI_IJK)
464
465
466 CALL FINL_ERR_MSG
467
468
469 RETURN
470 END SUBROUTINE SET_DEM_BCMI_IJK
471