File: N:\mfix\model\des\calc_collision_wall_mod.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CALC_COLLISION_WALL                                    C
4     !  Author: Rahul Garg                               Date: 1-Dec-2013   C
5     !                                                                      C
6     !  Purpose: subroutines for particle-wall collisions when cutcell is   C
7     !           used. Also contains rehack of routines for cfslide and     C
8     !           cffctow which might be different from the stand alone      C
9     !           routines. Eventually all the DEM routines will be          C
10     !           consolidated.                                              C
11     !                                                                      C
12     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
13           MODULE CALC_COLLISION_WALL
14     
15           PRIVATE
16           PUBLIC :: CALC_DEM_FORCE_WITH_WALL_STL,&
17           & CALC_DEM_THERMO_WITH_WALL_STL
18     
19           CONTAINS
20     
21     !VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV!
22     !                                                                      !
23     !                                                                      !
24     !                                                                      !
25     !                                                                      !
26     !                                                                      !
27     !                                                                      !
28     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
29           SUBROUTINE CALC_DEM_FORCE_WITH_WALL_STL
30     
31           USE run
32           USE param1
33           USE desgrid
34           USE discretelement
35           USE geometry
36           USE compar
37           USE constant
38           USE indices
39           USE stl
40           USE stl_functions_des
41           USE functions
42           Implicit none
43     
44           INTEGER :: LL
45           INTEGER :: IJK, NF
46           DOUBLE PRECISION ::OVERLAP_N, SQRT_OVERLAP
47     
48           DOUBLE PRECISION :: V_REL_TRANS_NORM, DISTSQ, RADSQ, CLOSEST_PT(DIMN)
49     ! local normal and tangential forces
50           DOUBLE PRECISION :: NORMAL(DIMN), VREL_T(DIMN), DIST(DIMN), DISTMOD
51           DOUBLE PRECISION, DIMENSION(DIMN) :: FTAN, FNORM, OVERLAP_T
52     
53           LOGICAL :: DES_LOC_DEBUG
54           INTEGER :: CELL_ID, cell_count
55           INTEGER :: PHASELL
56     
57           DOUBLE PRECISION :: TANGENT(DIMN)
58           DOUBLE PRECISION :: FTMD, FNMD
59     ! local values used spring constants and damping coefficients
60           DOUBLE PRECISION ETAN_DES_W, ETAT_DES_W, KN_DES_W, KT_DES_W
61     
62           double precision :: MAG_OVERLAP_T
63     
64           double precision :: line_t
65     ! flag to tell if the orthogonal projection of sphere center to
66     ! extended plane detects an overlap
67     
68           DOUBLE PRECISION :: MAX_DISTSQ, DISTAPART, FORCE_COH, R_LM
69           INTEGER :: MAX_NF, axis
70           DOUBLE PRECISION, DIMENSION(3) :: PARTICLE_MIN, PARTICLE_MAX, POS_TMP
71     
72           DES_LOC_DEBUG = .false. ;      DEBUG_DES = .false.
73           FOCUS_PARTICLE = -1
74     
75     !$omp parallel default(none) private(LL,ijk,MAG_OVERLAP_T,             &
76     !$omp    cell_id,radsq,particle_max,particle_min,tangent,              &
77     !$omp    axis,nf,closest_pt,dist,r_lm,distapart,force_coh,distsq,      &
78     !$omp    line_t,max_distsq,max_nf,normal,distmod,overlap_n,VREL_T,     &
79     !$omp    v_rel_trans_norm,phaseLL,sqrt_overlap,kn_des_w,kt_des_w,      &
80     !$omp    etan_des_w,etat_des_w,fnorm,overlap_t,ftan,ftmd,fnmd,pos_tmp) &
81     !$omp shared(max_pip,focus_particle,debug_des,                         &
82     !$omp    pijk,dg_pijk,i_of,j_of,k_of,des_pos_new,    &
83     !$omp    des_radius,facets_at_dg,vertex,  &
84     !$omp    hert_kwn,hert_kwt,kn_w,kt_w,des_coll_model_enum,mew_w,tow,    &
85     !$omp    des_etan_wall,des_etat_wall,dtsolid,fc,norm_face,             &
86     !$omp    wall_collision_facet_id,wall_collision_PFT,use_cohesion,      &
87     !$omp    van_der_waals,wall_hamaker_constant,wall_vdw_outer_cutoff,    &
88     !$omp    wall_vdw_inner_cutoff,asperities,surface_energy)
89     !$omp do
90           DO LL = 1, MAX_PIP
91     
92              IF(LL.EQ.FOCUS_PARTICLE) DEBUG_DES = .TRUE.
93     
94     ! skipping non-existent particles or ghost particles
95     ! make sure the particle is not classified as a new 'entering' particle
96     ! or is already marked as a potential exiting particle
97              IF(.NOT.IS_NORMAL(LL)) CYCLE
98     
99              CELL_ID = DG_PIJK(LL)
100     
101     ! If no neighboring facet in the surrounding 27 cells, then exit
102              IF(facets_at_dg(CELL_ID)%COUNT < 1) THEN
103                 WALL_COLLISION_FACET_ID(:,LL) = -1
104                 WALL_COLLISION_PFT(:,:,LL) = 0.0d0
105                 CYCLE
106              ENDIF
107     
108     ! Check particle LL for wall contacts
109              RADSQ = DES_RADIUS(LL)*DES_RADIUS(LL)
110     
111              particle_max(:) = des_pos_new( LL,:) + des_radius(LL)
112              particle_min(:) = des_pos_new( LL,:) - des_radius(LL)
113     
114              DO CELL_COUNT = 1, facets_at_dg(cell_id)%count
115     
116                 axis = facets_at_dg(cell_id)%dir(cell_count)
117     
118                 NF = facets_at_dg(cell_id)%id(cell_count)
119     
120     ! Compute particle-particle VDW cohesive short-range forces
121                 IF(USE_COHESION .AND. VAN_DER_WAALS) THEN
122     
123                    CALL ClosestPtPointTriangle(DES_POS_NEW(LL,:),          &
124                       VERTEX(:,:,NF), CLOSEST_PT(:))
125     
126                    DIST(:) = CLOSEST_PT(:) - DES_POS_NEW(LL,:)
127                    DISTSQ = DOT_PRODUCT(DIST, DIST)
128                    R_LM = 2*DES_RADIUS(LL)
129     
130                    IF(DISTSQ < (R_LM+WALL_VDW_OUTER_CUTOFF)**2) THEN
131                       IF(DISTSQ > (WALL_VDW_INNER_CUTOFF+R_LM)**2) THEN
132                          DistApart = (SQRT(DISTSQ)-R_LM)
133                          FORCE_COH = WALL_HAMAKER_CONSTANT*DES_RADIUS(LL) /&
134                             (12d0*DistApart**2)*(Asperities/(Asperities +  &
135                             DES_RADIUS(LL)) + ONE/(ONE+Asperities/         &
136                             DistApart)**2)
137                       ELSE
138                          FORCE_COH = 2d0*PI*SURFACE_ENERGY*DES_RADIUS(LL)* &
139                             (Asperities/(Asperities+DES_RADIUS(LL)) + ONE/ &
140                             (ONE+Asperities/WALL_VDW_INNER_CUTOFF)**2 )
141                       ENDIF
142                       FC(LL,:) = FC(LL,:) + DIST(:)*FORCE_COH/SQRT(DISTSQ)
143                    ENDIF
144                 ENDIF
145     
146                 if (facets_at_dg(cell_id)%min(cell_count) >    &
147                    particle_max(axis)) then
148                    call remove_collision(LL, nf, wall_collision_facet_id)
149                    cycle
150                 endif
151     
152                 if (facets_at_dg(cell_id)%max(cell_count) <    &
153                    particle_min(axis)) then
154                    call remove_collision(LL, nf, wall_collision_facet_id)
155                    cycle
156                 endif
157     
158     ! Checking all the facets is time consuming due to the expensive
159     ! separating axis test. Remove this facet from contention based on
160     ! a simple orthogonal projection test.
161     
162     ! Parametrize a line as p = p_0 + t normal and intersect with the
163     ! triangular plane. If t>0, then point is on the non-fluid side of
164     ! the plane. If the plane normal is assumed to point toward the fluid.
165     
166     ! -undefined, because non zero values will imply the sphere center
167     ! is on the non-fluid side of the plane. Since the testing
168     ! is with extended plane, this could very well happen even
169     ! when the particle is well inside the domain (assuming the plane
170     ! normal points toward the fluid). See the pic below. So check
171     ! only when line_t is negative
172     
173     !                 \   Solid  /
174     !                  \  Side  /
175     !                   \      /
176     !                    \    /
177     ! Wall 1, fluid side  \  /  Wall 2, fluid side
178     !                      \/
179     !                        o particle
180     !
181     ! line_t will be positive for wall 1 (incorrectly indicating center
182     ! is outside the domain) and line_t will be negative for wall 2.
183     !
184     ! Therefore, only stick with this test when line_t is negative and let
185     ! the separating axis test take care of the other cases.
186     
187     ! Since this is for checking static config, line's direction is the
188     ! same as plane's normal. For moving particles, the line's normal will
189     ! be along the point joining new and old positions.
190     
191                 line_t = DOT_PRODUCT(VERTEX(1,:,NF) - des_pos_new(LL,:),&
192                    NORM_FACE(:,NF))
193     
194     ! k - rad >= tol_orth, where k = -line_t, then orthogonal
195     ! projection is false. Substituting for k
196     ! => line_t + rad <= -tol_orth
197     ! choosing tol_orth = 0.01% of des_radius = 0.0001*des_radius
198     
199     ! Orthogonal projection will detect false positives even
200     ! when the particle does not overlap the triangle.
201     ! However, if the orthogonal projection shows no overlap, then
202     ! that is a big fat negative and overlaps are not possible.
203                 if((line_t.le.-1.0001d0*des_radius(LL))) then  ! no overlap
204                    call remove_collision(LL,nf,wall_collision_facet_id)
205                    CYCLE
206                 ENDIF
207     
208                 pos_tmp = DES_POS_NEW(LL,:)
209                 CALL ClosestPtPointTriangle(pos_tmp,             &
210                    VERTEX(:,:,NF), CLOSEST_PT(:))
211     
212                 DIST(:) = CLOSEST_PT(:) - DES_POS_NEW(LL,:)
213                 DISTSQ = DOT_PRODUCT(DIST, DIST)
214     
215                 IF(DISTSQ .GE. RADSQ - SMALL_NUMBER) THEN !No overlap exists
216                    call remove_collision(LL,nf,wall_collision_facet_id)
217                    CYCLE
218                 ENDIF
219     
220                 MAX_DISTSQ = DISTSQ
221                 MAX_NF = NF
222     
223     ! Assign the collision normal based on the facet with the
224     ! largest overlap.
225                 NORMAL(:) = DIST(:)/sqrt(DISTSQ)
226     
227     ! Facet's normal is correct normal only when the intersection is with
228     ! the face. When the intersection is with edge or vertex, then the
229     ! normal is based on closest pt and sphere center. The definition above
230     ! of the normal is generic enough to account for differences between
231     ! vertex, edge, and facet.
232     
233     ! Calculate the particle/wall overlap.
234                 DISTMOD = SQRT(MAX_DISTSQ)
235                 OVERLAP_N = DES_RADIUS(LL) - DISTMOD
236     
237     ! Calculate the translational relative velocity
238                 CALL CFRELVEL_WALL(LL, V_REL_TRANS_NORM,VREL_T,           &
239                    NORMAL, DISTMOD)
240     
241     ! Calculate the spring model parameters.
242                 phaseLL = PIJK(LL,5)
243     
244     ! Hertz vs linear spring-dashpot contact model
245                 IF (DES_COLL_MODEL_ENUM .EQ. HERTZIAN) THEN
246                    sqrt_overlap = SQRT(OVERLAP_N)
247                    KN_DES_W = hert_kwn(phaseLL)*sqrt_overlap
248                    KT_DES_W = hert_kwt(phaseLL)*sqrt_overlap
249                    sqrt_overlap = SQRT(sqrt_overlap)
250                    ETAN_DES_W = DES_ETAN_WALL(phaseLL)*sqrt_overlap
251                    ETAT_DES_W = DES_ETAT_WALL(phaseLL)*sqrt_overlap
252                 ELSE
253                    KN_DES_W = KN_W
254                    KT_DES_W = KT_W
255                    ETAN_DES_W = DES_ETAN_WALL(phaseLL)
256                    ETAT_DES_W = DES_ETAT_WALL(phaseLL)
257                 ENDIF
258     
259     ! Calculate the normal contact force
260                 FNORM(:) = -(KN_DES_W * OVERLAP_N * NORMAL(:) + &
261                    ETAN_DES_W * V_REL_TRANS_NORM * NORMAL(:))
262     
263     ! Calculate the tangential displacement.
264                 OVERLAP_T(:) = DTSOLID*VREL_T(:) + GET_COLLISION(LL,       &
265                    NF, WALL_COLLISION_FACET_ID, WALL_COLLISION_PFT)
266                 MAG_OVERLAP_T = sqrt(DOT_PRODUCT(OVERLAP_T, OVERLAP_T))
267     
268     ! Check for Coulombs friction law and limit the maximum value of the
269     ! tangential force on a particle in contact with a wall.
270                 IF(MAG_OVERLAP_T > 0.0) THEN
271     ! Tangential froce from spring.
272                    FTMD = KT_DES_W*MAG_OVERLAP_T
273     ! Max force before the on set of frictional slip.
274                    FNMD = MEW_W*sqrt(DOT_PRODUCT(FNORM,FNORM))
275     ! Direction of tangential force.
276                    TANGENT = OVERLAP_T/MAG_OVERLAP_T
277                    IF(FTMD < FNMD) THEN
278                       FTAN = -FTMD * TANGENT
279                    ELSE
280                       FTAN = -FNMD * TANGENT
281                       OVERLAP_T = (FNMD/KT_DES_W) * TANGENT
282                    ENDIF
283                 ELSE
284                    FTAN = 0.0
285                 ENDIF
286     ! Add in the tangential dashpot damping force
287                 FTAN = FTAN - ETAT_DES_W*VREL_T(:)
288     
289     ! Save the tangential displacement.
290                 CALL UPDATE_COLLISION(OVERLAP_T, LL, NF,                   &
291                    WALL_COLLISION_FACET_ID, WALL_COLLISION_PFT)
292     
293     ! Add the collision force to the total forces acting on the particle.
294                 FC(LL,:) = FC(LL,:) + FNORM(:) + FTAN(:)
295     
296     ! Add the torque force to the toal torque acting on the particle.
297                 TOW(LL,:) = TOW(LL,:) + DISTMOD*DES_CROSSPRDCT(NORMAL,FTAN)
298     
299              ENDDO
300     
301           ENDDO
302     !$omp end do
303     !$omp end parallel
304     
305           RETURN
306     
307            contains
308     
309     !......................................................................!
310     !  Function: GET_COLLISION                                             !
311     !                                                                      !
312     !  Purpose: Return the integrated (t0->t) tangential displacement.     !
313     !......................................................................!
314           FUNCTION GET_COLLISION(LLL,FACET_ID,WALL_COLLISION_FACET_ID,     &
315               WALL_COLLISION_PFT)
316     
317           use stl_dbg_des, only: write_this_stl
318           use stl_dbg_des, only: write_stls_this_dg
319     
320           use error_manager
321     
322           IMPLICIT NONE
323     
324           DOUBLE PRECISION :: GET_COLLISION(DIMN)
325           INTEGER, INTENT(IN) :: LLL,FACET_ID
326           INTEGER, INTENT(INOUT) :: WALL_COLLISION_FACET_ID(:,:)
327           DOUBLE PRECISION, INTENT(INOUT) :: WALL_COLLISION_PFT(:,:,:)
328           INTEGER :: CC, FREE_INDEX, LC, dgIJK
329     
330           INTEGER :: lSIZE1, lSIZE2, lSIZE3
331           INTEGER, ALLOCATABLE :: tmpI2(:,:)
332           DOUBLE PRECISION, ALLOCATABLE :: tmpR3(:,:,:)
333     
334           free_index = -1
335     
336           do cc = 1, COLLISION_ARRAY_MAX
337              if (facet_id == wall_collision_facet_id(cc,LLL)) then
338                 get_collision(:) = wall_collision_PFT(:,cc,LLL)
339                 return
340              else if (-1 == wall_collision_facet_id(cc,LLL)) then
341                 free_index = cc
342              endif
343           enddo
344     
345     ! Overwrite old data. This is needed because a particle moving from
346     ! one dg cell to another may no longer 'see' an STL before it moved
347     ! out of contact range. Therefore, the 'remove_collision' function
348     ! does not get called to cleanup the stale data.
349           if(-1 == free_index) then
350              dgIJK=DG_PIJK(LLL)
351              cc_lp: do cc=1, COLLISION_ARRAY_MAX
352                 do lc=1, facets_at_dg(dgIJK)%count
353                    if(wall_collision_facet_id(cc,LLL) == &
354                       facets_at_dg(dgIJK)%id(LC))  cycle cc_lp
355                 enddo
356                 free_index = cc
357                 exit cc_lp
358              enddo cc_lp
359           endif
360     
361     ! Last resort... grow the collision array
362           if(-1 == free_index) then
363              free_index=COLLISION_ARRAY_MAX+1
364              COLLISION_ARRAY_MAX = 2*COLLISION_ARRAY_MAX
365              CALL GROW_WALL_COLLISION(COLLISION_ARRAY_MAX)
366           endif
367     
368           wall_collision_facet_id(free_index,LLL) = facet_id
369           wall_collision_PFT(:,free_index,LLL) = ZERO
370           get_collision(:) = wall_collision_PFT(:,free_index,LLL)
371           return
372     
373           END FUNCTION GET_COLLISION
374     
375     
376     !......................................................................!
377     !  Subroutine: GROW_WALL_COLLISION                                     !
378     !                                                                      !
379     !  Purpose: Return the integrated (t0->t) tangential displacement.     !
380     !......................................................................!
381           SUBROUTINE GROW_WALL_COLLISION(NEW_SIZE)
382     
383           use discretelement
384     
385           use stl_dbg_des, only: write_this_stl
386           use stl_dbg_des, only: write_stls_this_dg
387     
388           use error_manager
389     
390           IMPLICIT NONE
391     
392           INTEGER, INTENT(IN) :: NEW_SIZE
393           INTEGER :: lSIZE1, lSIZE2, lSIZE3
394           INTEGER, ALLOCATABLE :: tmpI2(:,:)
395           DOUBLE PRECISION, ALLOCATABLE :: tmpR3(:,:,:)
396     
397           lSIZE1 = size(wall_collision_facet_id,1)
398           lSIZE2 = size(wall_collision_facet_id,2)
399     
400           allocate(tmpI2(NEW_SIZE, lSIZE2))
401           tmpI2(1:lSIZE1,:) = WALL_COLLISION_FACET_ID(1:lSIZE1,:)
402           call move_alloc(tmpI2, WALL_COLLISION_FACET_ID)
403           WALL_COLLISION_FACET_ID(lSIZE1+1:NEW_SIZE,:) = -1
404     
405           lSIZE1 = size(wall_collision_pft,1)
406           lSIZE2 = size(wall_collision_pft,2)
407           lSIZE3 = size(wall_collision_pft,3)
408     
409           allocate(tmpR3(lSIZE1, NEW_SIZE, lSIZE3))
410           tmpR3(:,1:lSIZE2,:) = WALL_COLLISION_PFT(:,1:lSIZE2,:)
411           call move_alloc(tmpR3, WALL_COLLISION_PFT)
412     
413           RETURN
414           END SUBROUTINE GROW_WALL_COLLISION
415     
416     
417     
418     
419     !......................................................................!
420     !  Function: UPDATE_COLLISION                                          !
421     !                                                                      !
422     !  Purpose: Update the integrated (t0->t) tangential displacement.     !
423     !......................................................................!
424           SUBROUTINE UPDATE_COLLISION(PFT, LLL, FACET_ID,                  &
425              WALL_COLLISION_FACET_ID, WALL_COLLISION_PFT)
426     
427           use error_manager
428           implicit none
429     
430           DOUBLE PRECISION, INTENT(IN) :: PFT(DIMN)
431           INTEGER, INTENT(IN) :: LLL,FACET_ID
432           INTEGER, INTENT(IN) :: WALL_COLLISION_FACET_ID(:,:)
433           DOUBLE PRECISION, INTENT(INOUT) :: WALL_COLLISION_PFT(:,:,:)
434           INTEGER :: CC
435     
436           do cc = 1, COLLISION_ARRAY_MAX
437              if (facet_id == wall_collision_facet_id(cc,LLL)) then
438                 wall_collision_PFT(:,cc,LLL) = PFT(:)
439                 return
440              endif
441           enddo
442     
443           CALL INIT_ERR_MSG("CALC_COLLISION_WALL_MOD: UPDATE_COLLISION")
444           WRITE(ERR_MSG, 1100)
445           CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
446     
447      1100 FORMAT('Error: COLLISION_ARRAY_MAX too small. ')
448     
449           END SUBROUTINE UPDATE_COLLISION
450     
451     !......................................................................!
452     !  Function: REMOVE_COLLISION                                          !
453     !                                                                      !
454     !  Purpose: Clear the integrated (t0->t) tangential displacement once  !
455     !  the collision is over (contact ended).                              !
456     !......................................................................!
457           SUBROUTINE REMOVE_COLLISION(LLL,FACET_ID,WALL_COLLISION_FACET_ID)
458     
459           use error_manager
460     
461           IMPLICIT NONE
462     
463           INTEGER, INTENT(IN) :: LLL,FACET_ID
464           INTEGER, INTENT(INOUT) :: WALL_COLLISION_FACET_ID(:,:)
465           INTEGER :: CC
466     
467           DO CC = 1, COLLISION_ARRAY_MAX
468              IF (FACET_ID == WALL_COLLISION_FACET_ID(CC,LLL)) THEN
469                 WALL_COLLISION_FACET_ID(CC,LLL) = -1
470                 RETURN
471              ENDIF
472           ENDDO
473     
474           END SUBROUTINE REMOVE_COLLISION
475     
476           END SUBROUTINE CALC_DEM_FORCE_WITH_WALL_STL
477     
478     
479     !----------------------------------------------------------------------!
480     !                                                                      !
481     !  Subroutine: CFRELVEL_WALL                                           !
482     !                                                                      !
483     !  Purpose: Calculate the normal and tangential components of the      !
484     !  relative velocity between a particle and wall contact.              !
485     !                                                                      !
486     !  Comments: Only the magnitude of the normal component is returned    !
487     !  whereas the full tangential vector is returned.                     !
488     !----------------------------------------------------------------------!
489           SUBROUTINE CFRELVEL_WALL(LL, VRN, VRT, NORM, DIST)
490     
491     ! Particle translational velocity
492           use discretelement, only: DES_VEL_NEW
493     ! Particle rotational velocity
494           use discretelement, only: OMEGA_NEW
495     ! Spatial array size (parameter)
496           use discretelement, only: DIMN
497     ! Function for calculating the cross prodcut
498           use discretelement, only: DES_CROSSPRDCT
499     
500           IMPLICIT NONE
501     
502     ! Dummy arguments:
503     !---------------------------------------------------------------------//
504     ! Particle index.
505           INTEGER, INTENT(IN) :: LL
506     ! Magnitude of the total relative translational velocity.
507           DOUBLE PRECISION, INTENT(OUT):: VRN
508     ! Total relative translational velocity (vector).
509           DOUBLE PRECISION, DIMENSION(DIMN), INTENT(OUT):: VRT
510     ! Unit normal from particle center to closest point on stl (wall)
511           DOUBLE PRECISION, DIMENSION(DIMN), INTENT(IN) :: NORM
512     ! Distance between particle center and stl (wall).
513           DOUBLE PRECISION, INTENT(IN) :: DIST
514     
515     ! Local variables
516     !---------------------------------------------------------------------//
517     ! Additional relative translational motion due to rotation
518           DOUBLE PRECISION, DIMENSION(DIMN) :: V_ROT
519     ! Total relative velocity at contact point
520           DOUBLE PRECISION, DIMENSION(DIMN) :: VRELTRANS
521     
522     ! Total relative velocity + rotational contribution
523           V_ROT = DIST*OMEGA_NEW(LL,:)
524           VRELTRANS(:) =  DES_VEL_NEW(LL,:) + DES_CROSSPRDCT(V_ROT, NORM)
525     
526     ! magnitude of normal component of relative velocity (scalar)
527           VRN = DOT_PRODUCT(VRELTRANS,NORM)
528     
529     ! total relative translational slip velocity at the contact point
530     ! Equation (8) in Tsuji et al. 1992
531           VRT(:) =  VRELTRANS(:) - VRN*NORM(:)
532     
533           RETURN
534           END SUBROUTINE CFRELVEL_WALL
535     
536     
537     !----------------------------------------------------------------//
538     ! SUBROUTINE: CALC_DEM_THERMO_WITH_WALL_STL
539     ! By: Aaron M.
540     ! Purpose: Compute heat transfer to particles due to boundaries
541     !----------------------------------------------------------------//
542     
543           SUBROUTINE CALC_DEM_THERMO_WITH_WALL_STL
544     
545           USE bc
546           USE compar, only: istart3, iend3, jstart3, jend3, kstart3, kend3
547           USE cutcell, only: cut_cell_at, area_cut, blocked_cell_at, small_cell_at
548           USE des_thermo
549           USE des_thermo_cond
550           USE desgrid
551           USE discretelement
552           USE exit, only: mfix_exit
553           USE functions
554           USE geometry, only: do_k, imax1, imax2, imin1, imin2
555           USE geometry, only: no_k, jmax1, jmax2, jmin1, jmin2
556           USE geometry, only: kmax1, kmax2, kmin1, kmin2, zlength
557           USE geometry, only: dx, dy, dz
558           USE param, only: dimension_i, dimension_j, dimension_k
559           USE param1
560           USE physprop
561           USE run, only: units, time
562           USE stl
563           USE stl_functions_des
564     
565           IMPLICIT NONE
566     
567           INTEGER :: LL             ! Loop index for particle
568           INTEGER :: CELL_ID        ! Desgrid cell index
569           INTEGER :: CELL_COUNT     ! Loop index for facets
570           INTEGER :: IJK_FLUID      ! IJK index for fluid cell
571           INTEGER :: I_FLUID, J_FLUID, K_FLUID
572           INTEGER :: I_Facet, J_Facet, K_Facet, IJK_Facet
573           INTEGER :: I1,J1,K1
574           INTEGER :: phase_LL       ! Phase index for particle
575     
576           INTEGER, PARAMETER :: MAX_CONTACTS = 12
577           INTEGER :: iFacet         ! loop index for facets
578           INTEGER :: count_Facets   ! counter for number of facets particle contacts
579           DOUBLE PRECISION, DIMENSION(3,MAX_CONTACTS) :: NORM_FAC_CONTACT
580           LOGICAL :: USE_FACET
581           LOGICAL :: DOMAIN_BDRY
582     
583           INTEGER :: NF             ! Facet ID
584           INTEGER :: AXIS           ! Facet direction
585           DOUBLE PRECISION :: RLENS_SQ ! lens radius squared
586           DOUBLE PRECISION :: RLENS ! lens radius
587           DOUBLE PRECISION :: RPART ! Particle radius
588     
589           ! vectors for minimum / maximum possible particle contacts
590           DOUBLE PRECISION, DIMENSION(3) :: PARTPOS_MIN, PARTPOS_MAX
591           DOUBLE PRECISION, DIMENSION(3) :: POS_TMP ! Temp vector
592           DOUBLE PRECISION, DIMENSION(3) :: DIST, NORMAL
593     
594           DOUBLE PRECISION, DIMENSION(3) :: CLOSEST_PT
595     
596           DOUBLE PRECISION :: line_t  ! Normal projection for part/wall
597           DOUBLE PRECISION :: DISTSQ ! Separation from particle and facet (squared)
598           DOUBLE PRECISION :: PROJ
599     
600           INTEGER :: BC_ID ! BC ID
601           INTEGER :: IBC   ! BC Loop Index
602     
603           DOUBLE PRECISION :: TWALL
604           DOUBLE PRECISION :: K_gas
605           DOUBLE PRECISION :: TPART
606           DOUBLE PRECISION :: OVERLAP
607           DOUBLE PRECISION :: QSWALL, AREA
608     
609           LOGICAL, SAVE :: OUTPUT_WARNING = .TRUE.
610     
611     
612     
613         !  IF(.NOT.DES_CONTINUUM_COUPLED.or.DES_EXPLICITLY_COUPLED)THEN
614         !     CALL PARTICLES_IN_CELL
615         !  ENDIF
616           DO LL = 1, MAX_PIP
617              ! Skip non-existent particles or ghost particles
618              IF (.NOT.IS_NORMAL(LL)) CYCLE
619              PHASE_LL = PIJK(LL,5)
620              IF(.NOT.CALC_COND_DES(PHASE_LL))CYCLE
621              ! Get desgrid cell index
622              CELL_ID = DG_PIJK(LL)
623     
624              ! Skip cells that do not have neighboring facet
625              IF (facets_at_dg(CELL_ID)%COUNT <1) CYCLE
626     
627              ! Store lens radius of particle
628              RLENS = (ONE+FLPC)*DES_RADIUS(LL)
629              RLENS_SQ = RLENS*RLENS
630     
631              RPART = DES_RADIUS(LL)
632     
633              ! Compute max/min particle locations
634              PARTPOS_MAX(:) = des_pos_new(LL,:) + RLENS
635              PARTPOS_MIN(:) = des_pos_new(LL,:) - RLENS
636     
637              ! Get fluid cell
638              I_FLUID = PIJK(LL,1)
639              J_FLUID = PIJK(LL,2)
640              K_FLUID = PIJK(LL,3)
641              IJK_FLUID = PIJK(LL,4)
642              TPART = DES_T_S(LL)
643     
644              ! Sometimes PIJK is not updated every solids timestep and
645              ! therefore doesn't give the correct fluid cell that
646              ! contains the particles.  If so, determine actual fluid
647              ! cell that contains the particle
648              IF(.NOT.DES_CONTINUUM_COUPLED.or.DES_EXPLICITLY_COUPLED)THEN
649                 IF(I_FLUID <= ISTART3 .OR. I_FLUID >= IEND3) THEN
650                    CALL PIC_SEARCH(I_FLUID, DES_POS_NEW(LL,1), XE,       &
651                    DIMENSION_I, IMIN2, IMAX2)
652                 ELSE
653                    IF((DES_POS_NEW(LL,1) >= XE(I_FLUID-1)) .AND.         &
654                    (DES_POS_NEW(LL,1) <  XE(I_FLUID))) THEN
655                    I_FLUID = I_FLUID
656                    ELSEIF((DES_POS_NEW(LL,1) >= XE(I_FLUID)) .AND.       &
657                       (DES_POS_NEW(LL,1) < XE(I_FLUID+1))) THEN
658                       I_FLUID = I_FLUID+1
659                    ELSEIF((DES_POS_NEW(LL,1) >= XE(I_FLUID-2)) .AND.     &
660                       (DES_POS_NEW(LL,1) < XE(I_FLUID-1))) THEN
661                       I_FLUID = I_FLUID-1
662                    ELSE
663                       CALL PIC_SEARCH(I_FLUID, DES_POS_NEW(LL,1), XE,    &
664                       DIMENSION_I, IMIN2, IMAX2)
665                    ENDIF
666                 ENDIF !(I)
667                 IF(J_FLUID <= JSTART3 .OR. J_FLUID >= JEND3) THEN
668                    CALL PIC_SEARCH(J_FLUID, DES_POS_NEW(LL,2), YN,       &
669                    DIMENSION_J, JMIN2, JMAX2)
670                 ELSE
671                    IF((DES_POS_NEW(LL,2) >= YN(J_FLUID-1)) .AND.         &
672                       (DES_POS_NEW(LL,2) <  YN(J_FLUID))) THEN
673                       J_FLUID = J_FLUID
674                    ELSEIF((DES_POS_NEW(LL,2) >= YN(J_FLUID)) .AND.       &
675                       (DES_POS_NEW(LL,2) < YN(J_FLUID+1))) THEN
676                       J_FLUID = J_FLUID+1
677                    ELSEIF((DES_POS_NEW(LL,2) >= YN(J_FLUID-2)) .AND.     &
678                       (DES_POS_NEW(LL,2) < YN(J_FLUID-1))) THEN
679                       J_FLUID = J_FLUID-1
680                    ELSE
681                       CALL PIC_SEARCH(J_FLUID, DES_POS_NEW(LL,2), YN,    &
682                       DIMENSION_J, JMIN2, JMAX2)
683                    ENDIF
684                 ENDIF !(J)
685     
686                 IF(NO_K) THEN
687                    K_FLUID = 1
688                 ELSE
689                    IF(K_FLUID <= KSTART3 .OR. K_FLUID >= KEND3) THEN
690                    CALL PIC_SEARCH(K_FLUID, DES_POS_NEW(LL,3), ZT,         &
691                       DIMENSION_K, KMIN2, KMAX2)
692                    ELSE
693                       IF((DES_POS_NEW(LL,3) >= ZT(K_FLUID-1)) .AND.        &
694                          (DES_POS_NEW(LL,3) < ZT(K_FLUID))) THEN
695                          K_FLUID = K_FLUID
696                       ELSEIF((DES_POS_NEW(LL,3) >= ZT(K_FLUID)) .AND.      &
697                          (DES_POS_NEW(LL,3) < ZT(K_FLUID+1))) THEN
698                          K_FLUID = K_FLUID+1
699                       ELSEIF((DES_POS_NEW(LL,3) >= ZT(K_FLUID-2)) .AND.    &
700                          (DES_POS_NEW(LL,3) < ZT(K_FLUID-1))) THEN
701                          K_FLUID = K_FLUID-1
702                       ELSE
703                          CALL PIC_SEARCH(K_FLUID, DES_POS_NEW(LL,3), ZT,   &
704                          DIMENSION_K, KMIN2, KMAX2)
705                       ENDIF
706                    ENDIF
707                 ENDIF !(K)
708                 IJK_FLUID = PIJK(LL,4)
709              ENDIF ! (NOT CONTINUUM_COUPLED OR EXPLICIT)
710     
711     
712     ! Initialize counter for number of facets
713              COUNT_FACETS = 0
714     
715              ! Loop over potential facets
716              DO CELL_COUNT = 1, facets_at_dg(CELL_ID)%count
717                 ! Get direction (axis) and facet id (NF)
718                 axis = facets_at_dg(cell_id)%dir(cell_count)
719                 NF = facets_at_dg(cell_id)%id(cell_count)
720     
721                 ! Check to see if facet is out of reach of particle
722                 if (facets_at_dg(cell_id)%min(cell_count) >    &
723                    partpos_max(axis)) then
724                    cycle
725                 endif
726                 if (facets_at_dg(cell_id)%max(cell_count) <    &
727                    partpos_min(axis)) then
728                    cycle
729                 endif
730     
731                 ! Compute projection (normal distance) from wall to particle
732                 line_t = DOT_PRODUCT(VERTEX(1,:,NF) - des_pos_new(LL,:),&
733                 &        NORM_FACE(:,NF))
734     
735                 ! If normal exceeds particle radius, particle is not in contact
736                 if((line_t.lt.-RLENS))CYCLE
737     
738                 ! Compute closest point on facet
739                 POS_TMP(:) = DES_POS_NEW(LL,:)
740                 CALL ClosestPtPointTriangle(POS_TMP, VERTEX(:,:,NF), &
741                 &    CLOSEST_PT(:))
742                 ! Compute position vector from particle to closest point on facet
743                 DIST(:) = CLOSEST_PT(:)-POS_TMP(:)
744                 DISTSQ = DOT_PRODUCT(DIST,DIST)
745     
746                 ! Skip particles that are more than lens radius from facet
747                 IF(DISTSQ .GE. (RLENS_SQ-SMALL_NUMBER))CYCLE
748     
749                 ! Only do heat transfer for particles that are directly above facet
750                 ! Heat transfer routines not generalized for edge contacts, but
751                 ! those should normally yield negligible heat transfer contributions
752     
753                 ! Normalize distance vector and compare to facet normal
754                 NORMAL(:)=-DIST(:)/sqrt(DISTSQ)
755                 PROJ = sqrt(abs(DOT_PRODUCT(NORMAL, NORM_FACE(:,NF))))
756                 IF(ABS(ONE-PROJ).gt.1.0D-6)CYCLE
757     
758                 ! Get overlap
759                 OVERLAP = RPART - SQRT(DISTSQ)
760     
761                 ! Initialize area for facet (for post-proc. flux)
762                 AREA = ZERO
763     
764                 ! Initialize BDRY FLAG
765                 DOMAIN_BDRY = .FALSE.
766                 ! Get BC_ID
767                 BC_ID = BC_ID_STL_FACE(NF)
768     
769                 ! BC_ID is set to 0 in pre-proc if stl is a domain boundary
770                 IF(BC_ID.eq.0)then
771                    I1=I_FLUID
772                    J1=J_FLUID
773                    K1=K_FLUID
774                    DOMAIN_BDRY = .TRUE.
775     
776                    ! Domain boundary, figure out real boundary ID
777                    IF(NORM_FACE(1,NF).ge.0.9999)THEN
778                       ! WEST face
779                       I1=IMIN2
780                       IF(DO_K)THEN
781                          AREA = DY(J_FLUID)*DZ(K_FLUID)
782                       ELSE
783                          AREA = DY(J_FLUID)*ZLENGTH
784                       ENDIF
785                    ELSEIF(NORM_FACE(1,NF).le.-0.9999)THEN
786                       ! EAST FACE
787                       I1=IMAX2
788                       IF(DO_K)THEN
789                          AREA = DY(J_FLUID)*DZ(K_FLUID)
790                       ELSE
791                          AREA = DY(J_FLUID)*ZLENGTH
792                       ENDIF
793                    ELSEIF(NORM_FACE(2,NF).ge.0.9999)THEN
794                       ! SOUTH FACE
795                       J1=JMIN2
796                       IF(DO_K)THEN
797                          AREA = DX(I_FLUID)*DZ(K_FLUID)
798                       ELSE
799                          AREA = DX(I_FLUID)*ZLENGTH
800                       ENDIF
801                    ELSEIF(NORM_FACE(2,NF).le.-0.9999)THEN
802                       ! NORTH FACE
803                       J1=JMAX2
804                       IF(DO_K)THEN
805                          AREA = DX(I_FLUID)*DZ(K_FLUID)
806                       ELSE
807                          AREA = DX(I_FLUID)*ZLENGTH
808                       ENDIF
809                    ELSEIF(NORM_FACE(3,NF).ge.0.9999)THEN
810                       ! BOTTOM FACE
811                       K1=KMIN2
812                       AREA = DX(I_FLUID)*DY(J_FLUID)
813     
814                    ELSEIF(NORM_FACE(3,NF).le.-0.9999)THEN
815                       ! TOP FACE
816                       K1=KMAX2
817                       AREA = DX(I_FLUID)*DY(J_FLUID)
818                    ELSE
819                       WRITE( *,*)'PROBLEM, COULD NOT FIND DOMAIN BOUNDARY'
820                       WRITE(*,*)' In calc_thermo_des_wall_stl'
821                       call mfix_exit(1)
822     
823                    ENDIF
824     
825     
826                    ! Loop through defined BCs to see which one particle neighbors
827                    DO IBC = 1, DIMENSION_BC
828                       IF(.NOT.BC_DEFINED(IBC))CYCLE
829                       IF (I1.ge.BC_I_W(IBC).and.I1.le.BC_I_E(IBC).and.&
830                           J1.ge.BC_J_S(IBC).and.J1.le.BC_J_N(IBC).and.&
831                           K1.ge.BC_K_B(IBC).and.K1.le.BC_K_T(IBC))THEN
832                           BC_ID = IBC
833                           exit
834                       ENDIF
835                    ENDDO
836                    IF(BC_ID.eq.0)then
837                       IF(OUTPUT_WARNING)THEN
838                          write(*,*)'WARNING: Could not find BC'
839                          write(*,*)'Check input deck to make sure domain boundaries &
840                          are defined'
841                          write(*,*)'SUPPRESSING FUTURE OUTPUT'
842                          write(*,*)'DES_POS_NEW'
843                          write(*,'(3(F12.5, 3X))')(DES_POS_NEW(LL,IBC),IBC=1,3)
844                          write(*,*)'I,J,K'
845                          write(*,*)I1,J1,K1
846                          write(*,*)'CLOSEST PT'
847                          write(*,'(3(F12.5, 3X))')(CLOSEST_PT(IBC),IBC=1,3)
848                          write(*,*)'NORM_FACE'
849                          write(*,'(3(F12.5, 3X))')(NORM_FACE(IBC,NF),IBC=1,3)
850                          OUTPUT_WARNING = .FALSE.
851                       ENDIF
852                       CYCLE  !
853                    ENDIF
854     
855                 ENDIF !Domain Boundary (facet ID was 0)
856     
857                 IF (BC_TYPE_ENUM(BC_ID) == NO_SLIP_WALL .OR. &
858                    BC_TYPE_ENUM(BC_ID) == FREE_SLIP_WALL .OR. &
859                    BC_TYPE_ENUM(BC_ID) == PAR_SLIP_WALL .OR. &
860                    BC_TYPE_ENUM(BC_ID) == CG_NSW .OR. &
861                    BC_TYPE_ENUM(BC_ID) == CG_FSW .OR. &
862                    BC_TYPE_ENUM(BC_ID) == CG_PSW) THEN
863     
864     
865     
866                    ! CHECK TO MAKE SURE FACET IS UNIQUE
867                    USE_FACET=.TRUE.
868                    DO IFACET=1,count_facets
869                       ! DO CHECK BY ENSURING NORMAL VECTOR IS NEARLY PARALLEL
870                       PROJ = sqrt(abs(DOT_PRODUCT(NORMAL, NORM_FAC_CONTACT(:,IFACET))))
871                       IF(ABS(ONE-PROJ).gt.1.0D-6)THEN
872                          USE_FACET=.FALSE.
873                          EXIT
874                       ENDIF
875                    ENDDO
876                    IF(.NOT.USE_FACET)CYCLE
877     
878                    ! FACET IS UNIQUE
879                    count_facets=count_facets+1
880                    NORM_FAC_CONTACT(:,count_facets)=NORMAL(:)
881     
882     ! Do heat transfer
883                    ! GET WALL TEMPERATURE
884                    TWALL = BC_TW_S(BC_ID,phase_LL)
885     
886                    ! GET GAS THERMAL CONDUCTIVITY
887                    if(k_g0.eq.UNDEFINED)then
888                       ! Compute gas conductivity as is done in calc_k_g
889                       ! But use average of particle and wall temperature to be gas temperature
890     
891                       K_Gas = 6.02D-5*SQRT(HALF*(TWALL+TPART)/300.D0) ! cal/(s.cm.K)
892                       ! 1 cal = 4.183925D0 J
893                       IF (UNITS == 'SI') K_Gas = 418.3925D0*K_Gas !J/s.m.K
894                    else
895                       K_Gas=k_g0
896                    endif
897                    IF(TWALL.eq.UNDEFINED)CYCLE
898                    QSWALL = DES_CONDUCTION_WALL(OVERLAP,K_s0(phase_LL), &
899                    &        K_s0(phase_LL),K_Gas,TWALL, TPART, RPART, &
900                    &        RLENS, phase_LL)
901     
902     
903                    Q_Source(LL) = Q_Source(LL)+QSWALL
904     
905                    ! BELOW CODE IS ONLY NECESSARY FOR OUTPUTING
906                    ! DATA.  Need to know fluid cell that contact
907                    ! point resides in so that wall flux can be
908                    ! output correctly.
909     
910                    I_FACET = I_FLUID
911                    J_FACET = J_FLUID
912                    K_FACET = K_FLUID
913     
914                    ! This checks to see if the contact point was NOT on a domain boundary
915     
916                    IF(.NOT.DOMAIN_BDRY)THEN
917                       IF(CLOSEST_PT(1) >= XE(I_FLUID))THEN
918                          I_FACET = MIN(IMAX1, I_FLUID+1)
919                       ELSEIF(CLOSEST_PT(1) < XE(I_FLUID-1))THEN
920                          I_FACET = MAX(IMIN1, I_FLUID-1)
921                       ENDIF
922     
923                       IF(CLOSEST_PT(2) >= YN(J_FLUID))THEN
924                          J_FACET = MIN(JMAX1, J_FLUID+1)
925                       ELSEIF(CLOSEST_PT(2) < YN(J_FLUID-1))THEN
926                          J_FACET = MAX(JMIN1, J_FLUID-1)
927                       ENDIF
928                       IF(DO_K)THEN
929                          IF(CLOSEST_PT(3) >= ZT(K_FLUID))THEN
930                             K_FACET = MIN(KMAX1, K_FLUID+1)
931                          ELSEIF(CLOSEST_PT(3) < ZT(K_FLUID-1))THEN
932                             K_FACET = MAX(KMIN1, K_FLUID-1)
933                          ENDIF
934                       ENDIF
935                       IJK_FACET=funijk(I_facet,J_facet,K_facet)
936                       AREA=AREA_CUT(IJK_FACET)
937     
938                    ! AREA is left undefined if contact was with cut-cell surface
939                       IF (AREA.eq.ZERO)then
940                          I_FACET = I_FLUID
941                          J_FACET = J_FLUID
942                          K_FACET = K_FLUID
943                          IJK_FACET=funijk(I_facet,J_facet,K_facet)
944                          IF(NORM_FACE(1,NF).ge.0.9999)THEN
945                          ! WEST face
946                             IF(DO_K)THEN
947                                AREA = DY(J_FACET)*DZ(K_FACET)
948                             ELSE
949                                AREA = DY(J_FACET)*ZLENGTH
950                             ENDIF
951                          ELSEIF(NORM_FACE(1,NF).le.-0.9999)THEN
952                          ! EAST FACE
953                             IF(DO_K)THEN
954                                AREA = DY(J_FACET)*DZ(K_FACET)
955                             ELSE
956                                AREA = DY(J_FACET)*ZLENGTH
957                             ENDIF
958                          ELSEIF(NORM_FACE(2,NF).ge.0.9999)THEN
959                          ! SOUTH FACE
960                             IF(DO_K)THEN
961                                AREA = DX(I_FACET)*DZ(K_FACET)
962                             ELSE
963                                AREA = DX(I_FACET)*ZLENGTH
964                             ENDIF
965                          ELSEIF(NORM_FACE(2,NF).le.-0.9999)THEN
966                          ! NORTH FACE
967                             IF(DO_K)THEN
968                                AREA = DX(I_FACET)*DZ(K_FACET)
969                             ELSE
970                                AREA = DX(I_FACET)*ZLENGTH
971                             ENDIF
972                          ELSEIF(NORM_FACE(3,NF).ge.0.9999)THEN
973                          ! BOTTOM FACE
974                             AREA = DX(I_FACET)*DY(J_FACET)
975                          ELSEIF(NORM_FACE(3,NF).le.-0.9999)THEN
976                          ! TOP FACE
977                             AREA = DX(I_FACET)*DY(J_FACET)
978                          ENDIF
979                       ENDIF ! Area==0 (because cut-cell facet exists in non cut-cell)
980                    ENDIF ! NOT DOMAIN_BDRY
981                    IJK_FACET = FUNIJK(I_FACET,J_FACET,K_FACET)
982                    ! AM error check
983                    IF(.NOT.FLUID_AT(IJK_FACET).AND. &
984                    &  .NOT.BLOCKED_CELL_AT(IJK_FACET))THEN
985                       write(*,*)'ERROR: Cell containing facet is not a fluid &
986                       &  cell or a blocked cell'
987                       write(*,*)FLUID_AT(IJK_FACET), BLOCKED_CELL_AT(IJK_FACET)
988                       write(*,*)'PART POS',(DES_POS_NEW(LL,IBC),IBC=1,3)
989                       write(*,*)'FACET NORM',(NORM_FACE(IBC,NF),IBC=1,3)
990                       write(*,*)'BC_ID', BC_ID
991                       write(*,*)'I,J,K (Facet)', I_FACET,J_FACET,K_FACET
992     
993                       call mfix_exit(1)
994                    ENDIF
995     
996                    IF(FLUID_AT(IJK_FACET))THEN
997                       DES_QW_Cond(IJK_FACET,phase_LL) = &
998                          DES_QW_Cond(IJK_FACET, phase_LL) + QSWALL/AREA
999                    ENDIF
1000     
1001                 ENDIF ! WALL BDRY
1002              ENDDO                  ! CELL_COUNT (facets)
1003           ENDDO  ! LL
1004           RETURN
1005     
1006           END SUBROUTINE CALC_DEM_THERMO_WITH_WALL_STL
1007     
1008           end module CALC_COLLISION_WALL
1009