File: /nfs/home/0/users/jenkins/mfix.git/model/des/set_bc_dem_mi.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: SET_BC_DEM_MI                                           !
4     !  Author: J.Musser                                   Date: 23-Nov-09  !
5     !                                                                      !
6     !  Purpose:                                                            !
7     !                                                                      !
8     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
9           SUBROUTINE SET_BC_DEM_MI
10     
11     !-----------------------------------------------
12     ! Modules
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     ! Local variables
36     !-----------------------------------------------
37           INTEGER BCV_I      ! BC loop counter
38           INTEGER M, MM           ! Mass phase loop counter
39           INTEGER HOLD, I         ! Dummy values
40           INTEGER RANGE_TOP, RANGE_BOT ! Dummy values
41           INTEGER PHASE_CNT        ! Number of solid phases at bc
42           INTEGER PHASE_LIST(DIM_M) ! List of phases used in current bc
43     
44     ! the number of particles injected in a solids time step
45           DOUBLE PRECISION NPMpSEC(DES_MMAX) ! For solid phase m
46           DOUBLE PRECISION NPpSEC
47           DOUBLE PRECISION NPpDT        ! Total for BC
48           DOUBLE PRECISION SCALED_VAL
49           DOUBLE PRECISION MAX_DIA ! Max diameter of incoming particles at bc
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     ! Temp inlet velocity for solids phase M
60           DOUBLE PRECISION VEL_TMP(DIM_M)
61           DOUBLE PRECISION EPs_TMP(DIM_M)
62     
63     ! Minimum/maximum solids velocity at inlet.  Also used in the iterative
64     ! steps as the starting and ending velocities
65           DOUBLE PRECISION  MIN_VEL, MAX_VEL
66     
67           DOUBLE PRECISION  MINIPV, MAXIPV
68     
69           INTEGER :: OCCUPANTS
70     ! jump_here
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     ! Allocate the MI data structures for NEW and RES2 cases. Allocation for
80     ! RES1 cases is done prior to reading the RES file.
81           IF(RUN_TYPE /= 'RESTART_1') CALL ALLOCATE_DEM_MI
82     
83     ! Loop over BCs that flagged for DEM mass inflow.
84           DO BCV_I = 1, DEM_BCMI
85     
86     ! Get the user defined BC ID.
87              BCV = DEM_BCMI_MAP(BCV_I)
88     
89              if(dFlag) write(*,"(2/,'Setting DEM_MI:',I3)") BCV_I
90     
91     ! The number of mass phases at this inlet.  While a system may be
92     ! polydisperse, the inlet could consist of a single mass phase
93              PHASE_CNT = 0
94     ! The mass phase indices of incoming particles at this inlet
95              PHASE_LIST(:) = -1
96     ! The max diameter of incoming particles at this inlet
97              MAX_DIA = ZERO
98     
99     ! Determine if the inlet is mono or polydisperse
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     ! Set the polydispersity flag.
109              DEM_MI(BCV_I)%POLYDISPERSE = (PHASE_CNT > 1)
110     
111     ! Layout the feed pattern.
112              CALL LAYOUT_MI_DEM(BCV, BCV_I, MAX_DIA)
113     
114     ! Initialize
115              MAX_VEL = ZERO
116              NPMpSEC(:) = ZERO
117     
118     ! Calculate the individual velocities for each solid phase
119              DO MM = 1, PHASE_CNT
120                 M = PHASE_LIST(MM)
121     
122     ! Pull off the BC velocity normal to the flow plane.
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     ! Check for min/max inlet velocity
130                 MAX_VEL = MAX(ABS(VEL_TMP(M)), MAX_VEL)
131     ! Calculate volumetric flow rate to convert to particle count. BC_AREA
132     ! was already corrected for cut cells and velocity was recalculated
133     ! to ensure user-specified mass or volumetric flow rates.
134                 VOL_FLOW = VEL_TMP(M) * BC_AREA(BCV) * BC_EP_S(BCV,M+SMAX)
135     ! Calculate the number of particles of mass phase M are injected per
136     ! second for each solid phase present at the boundary
137                 NPMpSEC(M) = VOL_FLOW / (PI/6.d0*DES_D_P0(M)**3)
138     ! Write some debugging information if needed.
139                 if(dFlag) write(*,1100) M, VEL_TMP(M), NPMpSEC(M)
140              ENDDO
141     
142     ! Total number of particles at BCV injected per second
143              NPpSEC = 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     ! The number of total particles per solid time step DTSOLID
152              NPpDT = NPpSEC * DTSOLID
153     
154     ! Inject one particle every PI_FACTOR solids time steps.
155              IF(NPpDT .LT. 1)THEN
156                 PI_COUNT(BCV_I) = 1
157                 PI_FACTOR(BCV_I) = FLOOR(real(1.d0/NPpDT))
158     ! Inject PI_COUNT particles every soilds time step.
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     ! Calculate the minimum inlet velocity. The cutoff is associated with
167     ! square packing of disks on a plane.
168              MINIPV = MAX_DIA/( DTSOLID*dble(PI_FACTOR(BCV_I)) * dble(     &
169                 FLOOR( real(OCCUPANTS)/real(PI_COUNT(BCV_I)))))
170     ! Calculate the velocity needed to ensure that half the inlet is free.
171     ! Inlets with velocities greater than this value can be randomly seeded,
172     ! otherwise, particles are seeded in according to the grid.
173              MAXIPV = 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     ! Set all BC solids velocities to the largest velocity and recalculate
195     ! BC_EP_s to determine the magnitude of the change.
196              EPs_ERR = 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     ! Over-write the current BC value.
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     ! If the net change in solids volume fraction is greatere than 0.01,
216     ! flag this as an error and exit. >> Let the user fix the input.
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     ! Report the amount of changes imposed on the BC in setting a
224     ! uniform inlet velocity.
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     ! For polydisperse inlets, construct the DES_POLY_LAYOUT array
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     ! For monodisperse inlets, store the single mass phase used
276              ELSE
277                 DEM_BC_POLY_LAYOUT(BCV_I,:) = PHASE_LIST(1)
278              ENDIF
279     
280     
281     ! Calculate des mass inlet time; time between injection.  If the run
282     ! type is RESTART_1, DES_MI_TIME will be picked up from the restart file
283     ! with an updated value.
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     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
312     !                                                                      !
313     !  Subroutine: SET_BC_DEM                                              !
314     !  Author: J.Musser                                   Date: 13-Jul-09  !
315     !                                                                      !
316     !  Purpose: Check the data provided for the des mass inflow boundary   !
317     !  condition and flag errors if the data is improper.  This module is  !
318     !  also used to convert the proveded information into the format       !
319     !  necessary for the dependent subrountines to function properly.      !
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     ! Loop over all inflow BCs to get an approximate count of the number
363     ! of fluid cells that are adjacent to them.
364           MAX_CELLS = 0
365           DO BCV_I=1, DEM_BCMI
366              BCV = DEM_BCMI_MAP(BCV_I)
367     
368     ! Set the search area a little bigger than the inlet area.
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     ! Allocate an array to hold the IJK values. This array should be
392     ! more than enough to store all the IJKs.
393           allocate( LOC_DEM_BCMI_IJK(MAX_CELLS) )
394     
395     
396     ! Loop over the IJKs for each BC and store only the IJKs that you
397     ! own as well as the start/end array positions for each BC.
398           LC = 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     ! Depending on the flow plane, the 'common' index needs set to reference
411     ! the fluid cell.
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     ! Store the IJKs.
429              DO K = K_b, K_t
430              DO J = J_s, J_n
431              DO I = I_w, I_e
432     ! Skip cells that this rank does not own or are considered dead.
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     ! Allocate the global store arrary array. This changes across MPI ranks.
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