File: RELATIVE:/../../../mfix.git/model/des/mass_inflow_dem.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: DES_MASS_INLET                                          !
4     !  Author: J.Musser                                   Date: 13-Jul-09  !
5     !                                                                      !
6     !  Purpose:  This routine fills in the necessary information for new   !
7     !  particles entereing the system.                                     !
8     !                                                                      !
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
10           SUBROUTINE MASS_INFLOW_DEM
11     
12           use bc
13           use des_allocate
14           use des_bc
15           use discretelement
16           use functions
17     
18           use mpi_utility, only: BCAST
19     
20           implicit none
21     
22           INTEGER :: IP, LS, M, NP, IJK, LC
23           INTEGER :: BCV, BCV_I
24           LOGICAL :: CHECK_FOR_ERRORS, OWNS
25     
26     ! I/J/K index of fluid cell containing the new particle.
27           INTEGER :: IJKP(3)
28     
29           DOUBLE PRECISION :: DIST, POS(3)
30     ! Random numbers -shared by all ranks-
31           DOUBLE PRECISION :: RAND(3)
32     
33           CHECK_FOR_ERRORS = .FALSE.
34     
35           DO BCV_I = 1, DEM_BCMI
36              BCV = DEM_BCMI_MAP(BCV_I)
37     
38              DO LC=DEM_BCMI_IJKSTART(BCV_I), DEM_BCMI_IJKEND(BCV_I)
39                IJK = DEM_BCMI_IJK(LC)
40     
41                 DO LS= 1,DG_PIC(IJK)%ISIZE
42                    NP = DG_PIC(IJK)%P(LS)
43                    IF(IS_EXITING(NP) .or. IS_EXITING_GHOST(NP)) CYCLE
44                    SELECT CASE (BC_PLANE(BCV))
45                    CASE('N'); DIST = DES_POS_NEW(2,NP) - YN(BC_J_s(BCV))
46                    CASE('S'); DIST = YN(BC_J_s(BCV)-1) - DES_POS_NEW(2,NP)
47                    CASE('E'); DIST = DES_POS_NEW(1,NP) - XE(BC_I_w(BCV))
48                    CASE('W'); DIST = XE(BC_I_w(BCV)-1) - DES_POS_NEW(1,NP)
49                    CASE('T'); DIST = DES_POS_NEW(3,NP) - ZT(BC_K_b(BCV))
50                    CASE('B'); DIST = ZT(BC_K_b(BCV)-1) - DES_POS_NEW(3,NP)
51                    END SELECT
52     ! The particle is still inside the domain
53                    IF(DIST > DES_RADIUS(NP)) CALL SET_NORMAL(NP)
54                 ENDDO
55              ENDDO
56     
57     ! All ranks generate random numbers but PE_IO BCASTS its values so that
58     ! the calculated position (with random wiggles) is the same on all ranks.
59              CALL RANDOM_NUMBER(RAND)
60              CALL BCAST(RAND)
61     
62     ! Check if any particles need seeded this time step.
63              IF(DEM_MI_TIME(BCV_I) > S_TIME) CYCLE
64     
65              LS = 1
66     
67     ! Loop over the particles being injected
68              PLoop: DO IP = 1, PI_COUNT(BCV_I)
69     
70     ! Increment the global max particle ID (all ranks).
71                 iMAX_GLOBAL_ID = iMAX_GLOBAL_ID + 1
72     
73     ! Determine the location and phase of the incoming particle.
74                 CALL SEED_NEW_PARTICLE(BCV, BCV_I, RAND, M, POS, IJKP, OWNS)
75     
76     ! Only the rank receiving the new particle needs to continue
77                 IF(.NOT.OWNS) CYCLE PLoop
78     
79     ! Increment the number of particle on the processor by one. If the max
80     ! number of particles is exceeded, set the error flag and cycle.
81                 PIP = PIP + 1
82                 CALL PARTICLE_GROW(PIP)
83     
84     ! Find the first free space in the particle existance array.
85                 NP_LP: DO NP = LS, MAX_PIP
86                    IF(IS_NONEXISTENT(NP)) THEN
87                       LS = NP
88                       EXIT NP_LP
89                    ENDIF
90                 ENDDO NP_LP
91     
92     ! Set the particle's global ID.
93                 iGLOBAL_ID(LS) = iMAX_GLOBAL_ID
94     
95     ! Set the properties of the new particle.
96                 CALL SET_NEW_PARTICLE_PROPS(BCV, M, LS, POS, IJKP)
97     
98              ENDDO PLoop
99     
100     ! Update the time for seeding the next particle.
101              DEM_MI_TIME(BCV_I) = S_TIME + PI_FACTOR(BCV_I)*DTSOLID
102     ! Set the flag for error checking.
103              CHECK_FOR_ERRORS = .TRUE.
104           ENDDO
105     
106           IF(CHECK_FOR_ERRORS) THEN
107           ENDIF
108     
109           RETURN
110           END SUBROUTINE MASS_INFLOW_DEM
111     
112     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
113     !                                                                      !
114     !  Subroutine: SEED_NEW_PARTICLE                                       !
115     !  Author: J.Musser                                   Date: 14-Aug-09  !
116     !                                                                      !
117     !  Purpose:  This routine uses the classification information to place !
118     !  a new particle in the proper location.                              !
119     !                                                                      !
120     !  Comments:                                                           !
121     !                                                                      !
122     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
123           SUBROUTINE SEED_NEW_PARTICLE(lBCV, lBCV_I, lRAND, lM, lPOS,      &
124              lIJKP, lOWNS)
125     
126     !-----------------------------------------------
127     ! Modules
128     !-----------------------------------------------
129           USE compar
130           USE bc
131           USE des_bc
132           USE desgrid
133           USE discretelement
134           USE funits
135           USE geometry
136           USE param1
137           USE physprop
138           IMPLICIT NONE
139     !-----------------------------------------------
140     ! Dummy arguments
141     !-----------------------------------------------
142     ! The associated bc no.
143           INTEGER, INTENT(IN) :: lBCV, lBCV_I
144     ! Random numbers
145           DOUBLE PRECISION, INTENT(IN) :: lRAND(3)
146     ! Phase of incoming particle.
147           INTEGER, INTENT(OUT) :: lM
148     ! Position of incoming particle.
149           DOUBLE PRECISION, INTENT(OUT) :: lPOS(3)
150     ! I/J/K index of fluid cell containing the new particle.
151           INTEGER, INTENT(OUT) :: lIJKP(3)
152     ! Logical indicating that the current rank is the owner
153           LOGICAL, INTENT(OUT) :: lOWNS
154     
155     !-----------------------------------------------
156     ! Local variables
157     !-----------------------------------------------
158     ! the associated bc no.
159     !      INTEGER :: BCV
160     ! Random position
161           DOUBLE PRECISION RAND1, RAND2
162     ! Array index of vacant position
163           INTEGER :: VACANCY
164     ! Number of array position.
165           INTEGER :: OCCUPANTS
166     
167           INTEGER :: RAND_I
168           INTEGER :: lI, lJ, lK
169     
170           DOUBLE PRECISION :: WINDOW
171     
172     !      IF(PARTICLE_PLCMNT(lBCV_I) == 'ORDR')THEN
173              VACANCY = DEM_MI(lBCV_I)%VACANCY
174              OCCUPANTS = DEM_MI(lBCV_I)%OCCUPANTS
175              DEM_MI(lBCV_I)%VACANCY = MOD(VACANCY,OCCUPANTS) + 1
176     !      ELSE
177     !         call bcast(lpar_rad)
178     !         call bcast(lpar_pos)
179     !         call bcast(m)
180     !      ENDIF
181     
182     
183     ! Assign the phase of the incoming particle.
184           IF(DEM_MI(lBCV_I)%POLYDISPERSE) THEN
185              RAND_I = ceiling(dble(NUMFRAC_LIMIT)*lRAND(1))
186              lM = DEM_BC_POLY_LAYOUT(lBCV_I, RAND_I)
187           ELSE
188              lM = DEM_BC_POLY_LAYOUT(lBCV_I,1)
189           ENDIF
190     
191           WINDOW = DEM_MI(lBCV_I)%WINDOW
192           RAND1 = HALF*DES_D_p0(lM) + (WINDOW - DES_D_p0(lM))*lRAND(1)
193           RAND2 = HALF*DES_D_p0(lM) + (WINDOW - DES_D_p0(lM))*lRAND(2)
194     
195     
196     ! Set the physical location and I/J/K location of the particle.
197           SELECT CASE(BC_PLANE(lBCV))
198           CASE('N','S')
199     
200              lPOS(1) = DEM_MI(lBCV_I)%P(VACANCY) + RAND1
201              lPOS(3) = DEM_MI(lBCV_I)%Q(VACANCY) + RAND2
202              lPOS(2) = DEM_MI(lBCV_I)%OFFSET
203     
204              lIJKP(1) = DEM_MI(lBCV_I)%W(VACANCY)
205              lIJKP(3) = DEM_MI(lBCV_I)%H(VACANCY)
206              lIJKP(2) = DEM_MI(lBCV_I)%L
207     
208           CASE('E','W')
209     
210              lPOS(2) = DEM_MI(lBCV_I)%P(VACANCY) + RAND1
211              lPOS(3) = DEM_MI(lBCV_I)%Q(VACANCY) + RAND2
212              lPOS(1) = DEM_MI(lBCV_I)%OFFSET
213     
214              lIJKP(2) = DEM_MI(lBCV_I)%W(VACANCY)
215              lIJKP(3) = DEM_MI(lBCV_I)%H(VACANCY)
216              lIJKP(1) = DEM_MI(lBCV_I)%L
217     
218           CASE('T','B')
219     
220              lPOS(1) = DEM_MI(lBCV_I)%P(VACANCY) + RAND1
221              lPOS(2) = DEM_MI(lBCV_I)%Q(VACANCY) + RAND2
222              lPOS(3) = DEM_MI(lBCV_I)%OFFSET
223     
224              lIJKP(1) = DEM_MI(lBCV_I)%W(VACANCY)
225              lIJKP(2) = DEM_MI(lBCV_I)%H(VACANCY)
226              lIJKP(3) = DEM_MI(lBCV_I)%L
227     
228           END SELECT
229     
230     ! Easier and cleaner to clear out the third component at the end.
231           IF(NO_K) lPOS(3) = ZERO
232     
233           lI = IofPOS(lPOS(1))
234           lJ = JofPOS(lPOS(2))
235           lK = KofPOS(lPOS(3))
236     
237           lOWNS = ((DG_ISTART <= lI) .AND. (lI <= DG_IEND) .AND.           &
238              (DG_JSTART <= lJ) .AND. (lJ <= DG_JEND))
239     
240           IF(DO_K) lOWNS = lOWNS .AND. (DG_KSTART<=lK) .AND. (lK<=DG_KEND)
241     
242           RETURN
243           END SUBROUTINE SEED_NEW_PARTICLE
244     
245     
246     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
247     !                                                                      !
248     !  Subroutine: SET_NEW_PARTICLE_PROPS                                  !
249     !  Author: J.Musser                                   Date: 14-Aug-09  !
250     !                                                                      !
251     !  Purpose:  Set the properties of the new particle.                   !
252     !                                                                      !
253     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
254           SUBROUTINE SET_NEW_PARTICLE_PROPS(lBCV, lM, lNP, lPOS, lIJKP)
255     
256           USE compar
257           USE bc
258           USE des_bc
259           USE discretelement
260           USE funits
261           USE geometry
262           USE param1
263           USE physprop
264           USE des_thermo
265           use des_rxns
266     
267           use run, only: ENERGY_EQ
268           use run, only: ANY_SPECIES_EQ
269           use constant, only: PI
270     
271           use desgrid
272           use indices
273           use functions
274     
275           IMPLICIT NONE
276     !-----------------------------------------------
277     ! Dummy arguments
278     !-----------------------------------------------
279     ! The associated bc no.
280           INTEGER, INTENT(IN) :: lBCV
281     ! Phase of incoming particle.
282           INTEGER, INTENT(IN) :: lM
283     ! Index of new particle
284           INTEGER, INTENT(IN) :: lNP
285     ! Position of incoming particle.
286           DOUBLE PRECISION, INTENT(IN) :: lPOS(3)
287     ! I/J/K index of fluid cell containing the new particle.
288           INTEGER, INTENT(IN) :: lIJKP(3)
289     ! I/J/K index of DES grid cell
290           INTEGER :: lI, lJ, lK
291     
292     ! Global phase index
293           INTEGER :: BC_M
294     
295     ! Shift the phase index by SMAX to refernece global variables.
296           BC_M = lM + SMAX
297     
298     ! The particle exists and is entering, not exiting nor a ghost particle
299           IF (IS_GHOST(lNP)) THEN
300              CALL SET_ENTERING_GHOST(lNP)
301           ELSE
302              CALL SET_ENTERING(lNP)
303           ENDIF
304     
305     ! Set the initial position values based on mass inlet class
306           DES_POS_NEW(:,lNP) = lPOS(:)
307           PPOS(:,lNP) = lPOS(:)
308           DES_VEL_NEW(1,lNP) = BC_U_s(lBCV,BC_M)
309           DES_VEL_NEW(2,lNP) = BC_V_s(lBCV,BC_M)
310           DES_VEL_NEW(3,lNP) = BC_W_s(lBCV,BC_M)
311     
312     ! Set the initial velocity values
313           IF (DO_OLD) THEN
314              DES_POS_OLD(:,lNP) = lPOS(:)
315              DES_VEL_OLD(:,lNP) = DES_VEL_NEW(:,lNP)
316              OMEGA_OLD(:,lNP) = 0
317           ENDIF
318     
319     ! Set the initial angular velocity values
320           OMEGA_NEW(:,lNP) = 0
321     
322     ! Set the initial angular position values
323           IF(PARTICLE_ORIENTATION) ORIENTATION(1:3,lNP) = INIT_ORIENTATION
324     
325     ! Set the particle radius value
326           DES_RADIUS(lNP) = HALF * DES_D_P0(lM)
327     
328     ! Set the particle density value
329           RO_Sol(lNP) = DES_RO_S(lM)
330     
331     ! Store the I/J/K indices of the particle.
332           PIJK(lNP,1:3) = lIJKP(:)
333           PIJK(lNP,4) = FUNIJK(lIJKP(1), lIJKP(2), lIJKP(3))
334     
335     ! Set the particle mass phase
336           PIJK(lNP,5) = lM
337     
338     ! Calculate the DES grid cell indices.
339           lI = min(DG_IEND2,max(DG_ISTART2,IOFPOS(DES_POS_NEW(1,lNP))))
340           lJ = min(DG_JEND2,max(DG_JSTART2,JOFPOS(DES_POS_NEW(2,lNP))))
341           IF(NO_K) THEN
342              lK = 1
343           ELSE
344              lK = min(DG_KEND2,max(DG_KSTART2,KOFPOS(DES_POS_NEW(3,lNP))))
345           ENDIF
346     ! Store the triple
347           DG_PIJK(lNP) = DG_FUNIJK(lI,lJ,lK)
348     
349     ! Calculate the new particle's Volume, Mass, OMOI
350           PVOL(lNP) = (4.0d0/3.0d0) * PI * DES_RADIUS(lNP)**3
351           PMASS(lNP) = PVOL(lNP) * RO_Sol(lNP)
352           OMOI(lNP) = 5.d0 / (2.d0 * PMASS(lNP) * DES_RADIUS(lNP)**2)
353     
354     ! Clear the drag force
355           IF(DES_EXPLICITLY_COUPLED) DRAG_FC(:,lNP) = ZERO
356     
357     ! If solving the energy equations, set the temperature
358           IF(ANY_SPECIES_EQ .OR. ENERGY_EQ ) THEN
359              DES_T_s_NEW(lNP) = BC_T_s(lBCV,BC_M)
360              DES_T_s_OLD(lNP) = DES_T_s_NEW(lNP)
361           ENDIF
362     
363     ! Set species mass fractions
364           IF((ENERGY_EQ .AND. C_PS0(BC_M)/=UNDEFINED) .OR. ANY_SPECIES_EQ)&
365              DES_X_s(lNP,1:NMAX(BC_M)) = BC_X_s(lBCV,BC_M,1:NMAX(BC_M))
366     
367     ! Calculate time dependent physical properties
368           CALL DES_PHYSICAL_PROP(lNP, .FALSE.)
369     
370     
371           RETURN
372           END SUBROUTINE SET_NEW_PARTICLE_PROPS
373     
374     
375     
376     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
377     !                                                                      !
378     !  Subroutine:  DES_NEW_PARTICLE_TEST                                  !
379     !                                                                      !
380     !  Purpose:  This routine checks if a new particle placed using the    !
381     !  random inlet was placed in contact with an existing particle.  If   !
382     !  so a flag is set indicating contact, and the new particle is        !
383     !  repositioned within the inlet domain.                               !
384     !                                                                      !
385     !  Author: J.Musser                                   Date: 14-Aug-09  !
386     !                                                                      !
387     !  Purpose: This routine has to be modified for parallel version       !
388     !           the parameter now accepts the lpar_rad and lpar_pos and    !
389     !           tests if it touches any particles                          !
390     !  Comments:                                                           !
391     !                                                                      !
392     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
393     
394           SUBROUTINE DES_NEW_PARTICLE_TEST(BCV_I,ppar_rad,ppar_pos,TOUCHING)
395     
396           USE compar
397           USE constant
398           USE des_bc
399           USE discretelement
400           USE funits
401           USE geometry
402           USE indices
403           USE param1
404           USE physprop
405           USE functions
406     
407           IMPLICIT NONE
408     !-----------------------------------------------
409     ! Dummy arguments
410     !-----------------------------------------------
411     ! index of boundary condition
412           INTEGER, INTENT(IN) :: BCV_I
413           DOUBLE PRECISION, INTENT(IN) :: ppar_pos(DIMN)
414           DOUBLE PRECISION, INTENT(IN) :: ppar_rad
415           LOGICAL, INTENT(INOUT) :: TOUCHING
416     !-----------------------------------------------
417     ! Local variables
418     !-----------------------------------------------
419     ! particle number id of a potential overlapping/contacting particle
420           INTEGER NP2
421     ! total number of particles in current ijk cell and loop counter
422           INTEGER NPG, LL
423     ! i, j, k indices along boundary used for loop counters
424           INTEGER I, J, K, IJK
425     ! for parallel processing
426           integer listart,liend,ljstart,ljend,lkstart,lkend
427     
428           DOUBLE PRECISION  DISTVEC(DIMN), DIST, R_LM
429     !-----------------------------------------------
430     
431           TOUCHING = .FALSE.
432     
433     ! For parallel processing the arrays has to be limited
434     !      select case (des_mi_class(bcv_i))
435     !      case ('XW','XE', 'YZw','YZe')
436     !         listart = gs_array(bcv_i,1)
437     !         liend = gs_array(bcv_i,2)
438     !         ljstart = max(gs_array(bcv_i,3),jstart)
439     !         ljend = min(gs_array(bcv_i,4),jend)
440     !         lkstart = max(gs_array(bcv_i,5),jstart)
441     !         lkend = min(gs_array(bcv_i,6),jend)
442     !      case ('YN','YS', 'XZn','XZs')
443     !         listart = max(gs_array(bcv_i,1),istart)
444     !         liend = min(gs_array(bcv_i,2),iend)
445     !         ljstart = gs_array(bcv_i,3)
446     !         ljend = gs_array(bcv_i,4)
447     !         lkstart = max(gs_array(bcv_i,5),jstart)
448     !         lkend = min(gs_array(bcv_i,6),jend)
449     !      case ('ZT','ZB', 'XYt','XYb')
450     !         listart = max(gs_array(bcv_i,1),istart)
451     !         liend = min(gs_array(bcv_i,2),iend)
452     !         ljstart = max(gs_array(bcv_i,3),jstart)
453     !         ljend = min(gs_array(bcv_i,4),jend)
454     !         lkstart = gs_array(bcv_i,5)
455     !         lkend = gs_array(bcv_i,6)
456     !      end select
457     
458           DO k = lkstart,lkend
459           DO j = ljstart,ljend
460           DO i = listart,liend
461     !      DO K = GS_ARRAY(BCV_I,5), GS_ARRAY(BCV_I,6)
462     !         DO J = GS_ARRAY(BCV_I,3), GS_ARRAY(BCV_I,4)
463     !           DO I =  GS_ARRAY(BCV_I,1), GS_ARRAY(BCV_I,2)
464                  IJK = FUNIJK(I,J,K)
465                  IF(ASSOCIATED(PIC(IJK)%P)) THEN
466                    NPG =  SIZE(PIC(IJK)%P)
467                    DO LL = 1, NPG
468                       NP2 = PIC(IJK)%P(LL)
469                       DISTVEC(:) = ppar_pos(:) - DES_POS_NEW(:,NP2)
470                       DIST = DOT_PRODUCT(DISTVEC,DISTVEC)
471                       R_LM = ppar_rad + DES_RADIUS(NP2)
472                       IF(DIST .LE. R_LM*R_LM) TOUCHING = .TRUE.
473                    ENDDO
474                  ENDIF
475                ENDDO
476              ENDDO
477            ENDDO
478     
479           RETURN
480           END SUBROUTINE DES_NEW_PARTICLE_TEST
481