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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CALC_COLLISION_WALL                                    C
4     !                                                                      C
5     !  Purpose: subroutines for particle-wall collisions when cutcell is   C
6     !           used. Also contains rehack of routines for cfslide and     C
7     !           cffctow which might be different from the stand alone      C
8     !           routines. Eventually all the DEM routines will be          C
9     !           consolidated.                                              C
10     !                                                                      C
11     !                                                                      C
12     ! CALC_FORCE_WITH_WALL_CUTFACE_STL: Particle-wall driver routine when  C
13     !    stl files are used                                                C
14     !                                                                      C
15     !  Author: Rahul Garg                               Date: 1-Dec-2013   C
16     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
17     
18           module CALC_COLLISION_WALL
19     
20           PRIVATE
21           PUBLIC:: CHECK_IF_PARTICLE_OVELAPS_STL, CALC_DEM_FORCE_WITH_WALL_STL, ADD_FACET
22     
23           CONTAINS
24     
25     
26     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
27     !                                                                      !
28     !  Subroutine: ADD_FACET                                               !
29     !                                                                      !
30     !                                                                      !
31     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
32           SUBROUTINE ADD_FACET(CELL_ID, FACET_ID)
33     
34           Use discretelement
35           USE stl
36     
37           implicit none
38     
39           INTEGER, INTENT(IN) :: cell_id, facet_id
40     
41           INTEGER, DIMENSION(:), ALLOCATABLE :: int_tmp
42           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: real_tmp
43     
44           INTEGER :: lSIZE2, ii
45           DOUBLE PRECISION :: smallest_extent, min_temp, max_temp
46     
47           IF(STL_FACET_TYPE(facet_id) /= FACET_TYPE_NORMAL) RETURN
48     
49           DO II = 1, CELLNEIGHBOR_FACET_NUM(CELL_ID)
50              IF(FACET_ID .EQ. CELLNEIGHBOR_FACET(CELL_ID)%P(II)) RETURN
51           ENDDO
52     
53           CELLNEIGHBOR_FACET_NUM(CELL_ID) = &
54              CELLNEIGHBOR_FACET_NUM(CELL_ID) + 1
55     
56           NO_NEIGHBORING_FACET_DES(CELL_ID)  = .FALSE.
57     
58           IF(cellneighbor_facet_num(cell_id) > &
59              cellneighbor_facet_max(cell_id)) THEN
60     
61              cellneighbor_facet_max(cell_id) = &
62              2*cellneighbor_facet_max(cell_id)
63     
64              lSIZE2 = size(cellneighbor_facet(cell_id)%p)
65              allocate(int_tmp(cellneighbor_facet_max(cell_id)))
66              int_tmp(1:lSIZE2) = cellneighbor_facet(cell_id)%p(1:lSIZE2)
67              call move_alloc(int_tmp,cellneighbor_facet(cell_id)%p)
68     
69              lSIZE2 = size(cellneighbor_facet(cell_id)%extentdir)
70              allocate(int_tmp(cellneighbor_facet_max(cell_id)))
71              int_tmp(1:lSIZE2) = &
72                 cellneighbor_facet(cell_id)%extentdir(1:lSIZE2)
73              call move_alloc(int_tmp,cellneighbor_facet(cell_id)%extentdir)
74     
75              lSIZE2 = size(cellneighbor_facet(cell_id)%extentmin)
76              allocate(real_tmp(cellneighbor_facet_max(cell_id)))
77              real_tmp(1:lSIZE2) = &
78                 cellneighbor_facet(cell_id)%extentmin(1:lSIZE2)
79              call move_alloc(real_tmp,cellneighbor_facet(cell_id)%extentmin)
80     
81              lSIZE2 = size(cellneighbor_facet(cell_id)%extentmax)
82              allocate(real_tmp(cellneighbor_facet_max(cell_id)))
83              real_tmp(1:lSIZE2) = &
84                 cellneighbor_facet(cell_id)%extentmax(1:lSIZE2)
85              call move_alloc(real_tmp,cellneighbor_facet(cell_id)%extentmax)
86     
87           ENDIF
88     
89           CELLNEIGHBOR_FACET(CELL_ID)%&
90              P(CELLNEIGHBOR_FACET_NUM(CELL_ID)) = FACET_ID
91           SMALLEST_EXTENT = HUGE(0.0)
92     
93           DO II=1,3
94              MIN_TEMP = MINVAL(VERTEX(:,II,FACET_ID))
95              MAX_TEMP = MAXVAL(VERTEX(:,II,FACET_ID))
96     
97              IF(ABS(MAX_TEMP - MIN_TEMP) < SMALLEST_EXTENT ) THEN
98                  CELLNEIGHBOR_FACET(CELL_ID)%&
99                     EXTENTDIR(CELLNEIGHBOR_FACET_NUM(CELL_ID)) = II
100                  CELLNEIGHBOR_FACET(CELL_ID)%&
101                     EXTENTMIN(CELLNEIGHBOR_FACET_NUM(CELL_ID)) = MIN_TEMP
102                  CELLNEIGHBOR_FACET(CELL_ID)%&
103                     EXTENTMAX(CELLNEIGHBOR_FACET_NUM(CELL_ID)) = MAX_TEMP
104                  SMALLEST_EXTENT = ABS(MAX_TEMP - MIN_TEMP)
105              ENDIF
106           ENDDO
107     
108           END SUBROUTINE ADD_FACET
109     
110     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
111     !                                                                      C
112     !  SUBROUTINE: CHECK_IF_PARTICLE_OVELAPS_STL                           C
113     !                                                                      C
114     !  Purpose: This subroutine is special written to check if a particle  C
115     !          overlaps any of the STL faces. The routine exits on         C
116     !          detecting an overlap. It is called after initial            C
117     !          generation of lattice configuration to remove out of domain C
118     !          particles                                                   C
119     !                                                                      C
120     !  Authors: Rahul Garg                               Date: 21-Mar-2014 C
121     !                                                                      C
122     !                                                                      C
123     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
124           SUBROUTINE CHECK_IF_PARTICLE_OVELAPS_STL(POSITION, RADIUS, &
125              OVERLAP_EXISTS)
126     
127           USE run
128           USE param1
129           USE discretelement, only: dimn, xe, yn, zt
130           USE geometry
131           USE constant
132           USE cutcell
133           USE indices
134           USE stl
135           USE compar
136           USE des_stl_functions
137           use desgrid
138           USE functions
139           Implicit none
140     
141           DOUBLE PRECISION, INTENT(IN) :: POSITION(DIMN), RADIUS
142           LOGICAL, INTENT(OUT) :: OVERLAP_EXISTS
143     
144           INTEGER I, J, K, IJK, NF
145     
146           DOUBLE PRECISION :: RADSQ, DISTSQ, DIST(DIMN), CLOSEST_PT(DIMN)
147           INTEGER :: COUNT_FAC, COUNT, contact_facet_count, NEIGH_CELLS, &
148           NEIGH_CELLS_NONNAT, &
149           LIST_OF_CELLS(27), CELL_ID, I_CELL, J_CELL, K_CELL, cell_count , &
150           IMINUS1, IPLUS1, JMINUS1, JPLUS1, KMINUS1, KPLUS1, PHASELL, LOC_MIN_PIP, &
151           LOC_MAX_PIP!, focus_particle
152     
153     
154           OVERLAP_EXISTS = .FALSE.
155     
156           I_CELL = MIN(DG_IEND2,MAX(DG_ISTART2,IOFPOS(POSITION(1))))
157           J_CELL = MIN(DG_JEND2,MAX(DG_JSTART2,JOFPOS(POSITION(2))))
158           K_CELL = 1
159           IF(DO_K) K_CELL =MIN(DG_KEND2,MAX(DG_KSTART2,KOFPOS(POSITION(3))))
160     
161           CELL_ID = DG_FUNIJK(I_CELL, J_CELL, K_CELL)
162           IF (NO_NEIGHBORING_FACET_DES(CELL_ID)) RETURN
163     
164           LIST_OF_CELLS(:) = -1
165           NEIGH_CELLS = 0
166           NEIGH_CELLS_NONNAT  = 0
167     
168           COUNT_FAC = LIST_FACET_AT_DES(CELL_ID)%COUNT_FACETS
169     
170           RADSQ = RADIUS*RADIUS
171     
172     ! Add the facets in the cell the particle currently resides in
173           IF (COUNT_FAC.gt.0)   then
174              NEIGH_CELLS = NEIGH_CELLS + 1
175              LIST_OF_CELLS(NEIGH_CELLS) = CELL_ID
176           ENDIF
177     
178           IPLUS1  =  MIN (I_CELL + 1, DG_IEND2)
179           IMINUS1 =  MAX (I_CELL - 1, DG_ISTART2)
180     
181           JPLUS1  =  MIN (J_CELL + 1, DG_JEND2)
182           JMINUS1 =  MAX (J_CELL - 1, DG_JSTART2)
183     
184           KPLUS1  =  MIN (K_CELL + 1, DG_KEND2)
185           KMINUS1 =  MAX (K_CELL - 1, DG_KSTART2)
186     
187           DO K = KMINUS1, KPLUS1
188              DO J = JMINUS1, JPLUS1
189                 DO I = IMINUS1, IPLUS1
190     
191                    IF(.NOT.dg_is_ON_myPE_plus1layers(I,J,K)) CYCLE
192     
193                    IJK = DG_FUNIJK(I,J,K)
194                    COUNT_FAC = LIST_FACET_AT_DES(IJK)%COUNT_FACETS
195     
196                    IF(COUNT_FAC.EQ.0) CYCLE
197                    distsq = zero
198     
199                    IF(POSITION(1) > XE(I)) DISTSQ = DISTSQ &
200                       + (POSITION(1)-XE(I))*(POSITION(1)-XE(I))
201     
202                    IF(POSITION(1) < XE(I) - DX(I)) DISTSQ = DISTSQ &
203                    + (XE(I) - DX(I) - POSITION(1))*(XE(I) - DX(I) - POSITION(1))
204     
205                    IF(POSITION(2) > YN(J)) DISTSQ = DISTSQ &
206                    + (POSITION(2)-YN(J))* (POSITION(2)-YN(J))
207     
208                    IF(POSITION(2) < YN(J) - DY(J)) DISTSQ = DISTSQ &
209                    + (YN(J) - DY(J) - POSITION(2))* (YN(J) - DY(J) - POSITION(2))
210     
211                    IF(POSITION( 3) > ZT(K)) DISTSQ = DISTSQ &
212                    + (POSITION(3)-ZT(K))*(POSITION(3)-ZT(K))
213     
214                    IF(POSITION(3) < ZT(K) - DZ(K)) DISTSQ = DISTSQ &
215                    + (ZT(K) - DZ(K) - POSITION(3))*(ZT(K) - DZ(K) - POSITION(3))
216                    IF (DISTSQ < RADSQ) then
217                       NEIGH_CELLS_NONNAT = NEIGH_CELLS_NONNAT + 1
218                       NEIGH_CELLS = NEIGH_CELLS + 1
219                       LIST_OF_CELLS(NEIGH_CELLS) = IJK
220                                     !WRITE(*,'(A10, 4(2x,i5))') 'WCELL  = ', IJK, I,J,K
221                    ENDIF
222                 ENDDO
223              ENDDO
224           ENDDO
225     
226           CONTACT_FACET_COUNT = 0
227     
228           DO CELL_COUNT = 1, NEIGH_CELLS
229              IJK = LIST_OF_CELLS(CELL_COUNT)
230     
231              DO COUNT = 1, LIST_FACET_AT_DES(IJK)%COUNT_FACETS
232                 NF = LIST_FACET_AT_DES(IJK)%FACET_LIST(COUNT)
233     
234                 CALL ClosestPtPointTriangle(POSITION(:), VERTEX(:,:,NF), CLOSEST_PT(:))
235     
236                 DIST(:) = POSITION(:) - CLOSEST_PT(:)
237                 DISTSQ = DOT_PRODUCT(DIST, DIST)
238     
239                 IF(DISTSQ .GE. RADSQ) CYCLE !No overlap exists, move on to the next facet
240     
241                 !Overlap detected
242                 !Set overlap_exists to true and exit
243                 OVERLAP_EXISTS = .true.
244                 RETURN
245     
246              ENDDO
247     
248           end DO
249     
250           RETURN
251     
252           END SUBROUTINE CHECK_IF_PARTICLE_OVELAPS_STL
253     
254     
255     
256     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
257     !                                                                      !
258     !                                                                      !
259     !                                                                      !
260     !                                                                      !
261     !                                                                      !
262     !                                                                      !
263     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
264           SUBROUTINE CALC_DEM_FORCE_WITH_WALL_STL
265     
266           USE run
267           USE param1
268           USE desgrid
269           USE discretelement
270           USE geometry
271           USE compar
272           USE constant
273           USE indices
274           USE stl
275           USE des_stl_functions
276           USE functions
277           Implicit none
278     
279           INTEGER :: LL
280           INTEGER IJK, NF
281           DOUBLE PRECISION OVERLAP_N, SQRT_OVERLAP
282     
283           DOUBLE PRECISION V_REL_TRANS_NORM, DISTSQ, RADSQ, CLOSEST_PT(DIMN)
284     ! local normal and tangential forces
285           DOUBLE PRECISION FNS1(DIMN), FNS2(DIMN)
286           DOUBLE PRECISION FTS1(DIMN), FTS2(DIMN)
287           DOUBLE PRECISION NORMAL(DIMN), TANGENT(DIMN), DIST(DIMN), DISTMOD
288           DOUBLE PRECISION, DIMENSION(DIMN) :: FTAN, FNORM, OVERLAP_T, PFT
289     
290           LOGICAL :: DES_LOC_DEBUG, PARTICLE_SLIDE
291           INTEGER :: COUNT_FAC, &
292           LIST_OF_CELLS(27), CELL_ID, I_CELL, J_CELL, K_CELL, cell_count
293           INTEGER :: IMINUS1, IPLUS1, JMINUS1, JPLUS1, KMINUS1, KPLUS1, PHASELL
294     
295           DOUBLE PRECISION :: CROSSP(DIMN)
296           DOUBLE PRECISION :: FTMD, FNMD
297     ! local values used spring constants and damping coefficients
298           DOUBLE PRECISION ETAN_DES_W, ETAT_DES_W, KN_DES_W, KT_DES_W
299     
300           double precision :: line_t
301           !flag to tell if the orthogonal projection of sphere center to
302           !extended plane detects an overlap
303           INTEGER, Parameter :: MAX_FACET_CONTS = 200
304     
305           DOUBLE PRECISION :: FORCE_HISTORY(DIMN), DTSOLID_TMP
306     
307     
308           DOUBLE PRECISION :: MAX_DISTSQ, DISTAPART, FORCE_COH, R_LM
309           INTEGER :: MAX_NF, axis
310           DOUBLE PRECISION, DIMENSION(3) :: PARTICLE_MIN, PARTICLE_MAX
311     
312           DES_LOC_DEBUG = .false. ;      DEBUG_DES = .false.
313           FOCUS_PARTICLE = -1
314     
315     !$omp parallel default(none) private(LL,ijk,count_fac,fts1,fts2,fns1,fns2,list_of_cells,  &
316     !$omp          cell_id,radsq,particle_max,particle_min,axis,nf,closest_pt,dist,r_lm,      &
317     !$omp          distapart,force_coh,distsq,line_t,max_distsq,max_nf,                       &
318     !$omp          normal,distmod,overlap_n,v_rel_trans_norm,tangent,phaseLL,sqrt_overlap,    &
319     !$omp          kn_des_w,kt_des_w,etan_des_w,etat_des_w,fnorm,overlap_t,                   &
320     !$omp          force_history,ftan,particle_slide,ftmd,fnmd,crossp,pft)              &
321     !$omp    shared(max_pip,focus_particle,debug_des,pea,no_neighboring_facet_des,pijk, &
322     !$omp           dg_pijk,list_facet_at_des,i_of,j_of,k_of,des_pos_new,des_radius,    &
323     !$omp           cellneighbor_facet_num,cellneighbor_facet,vertex,hert_kwn,hert_kwt, &
324     !$omp           kn_w,kt_w,des_coll_model_enum,mew_w,tow,   &
325     !$omp           des_etan_wall,des_etat_wall,dtsolid,dtsolid_tmp,fc,norm_face,wall_collision_facet_id,wall_collision_PFT,       &
326     !$omp           use_cohesion,van_der_waals,wall_hamaker_constant,wall_vdw_outer_cutoff,wall_vdw_inner_cutoff, &
327     !$omp           asperities,surface_energy)
328     !$omp do
329           DO LL = 1, MAX_PIP
330     
331              IF(LL.EQ.FOCUS_PARTICLE) DEBUG_DES = .TRUE.
332     
333     ! skipping non-existent particles or ghost particles
334              IF(.NOT.PEA(LL,1) .OR. PEA(LL,4)) CYCLE
335     
336     !---------------------------------------------------------------------
337     ! make sure the particle is not classified as a new 'entering' particle
338     ! or is already marked as a potential exiting particle
339     
340              IF( PEA(LL,2) .OR. PEA(LL,3)) CYCLE
341     
342     ! If no neighboring facet in the surrounding 27 cells, then exit
343              IF (NO_NEIGHBORING_FACET_DES(DG_PIJK(LL)))  cycle
344     
345             IF(DEBUG_DES.AND.LL.EQ.FOCUS_PARTICLE) THEN
346                 IJK = PIJK(LL,4)
347                 COUNT_FAC = LIST_FACET_AT_DES(IJK)%COUNT_FACETS
348     
349                 WRITE(*,*) 'NUMBER OF FACETS = ', I_OF(IJK), J_OF(IJK), K_OF(IJK), IJK
350                 WRITE(*,*) 'NUMBER OF FACETS = ', COUNT_FAC, I_OF(IJK), J_OF(IJK), K_OF(IJK)
351     
352                 WRITE(*,'(A, 3(2x, g17.8))') 'POS = ', DES_POS_NEW(:, LL)
353              ENDIF
354     
355              FTS1(:) = ZERO
356              FTS2(:) = ZERO
357              FNS1(:) = ZERO
358              FNS2(:) = ZERO
359     
360     ! Check particle LL for wall contacts
361     
362              LIST_OF_CELLS(:) = -1
363     
364              CELL_ID = DG_PIJK(LL)
365     
366              COUNT_FAC = LIST_FACET_AT_DES(CELL_ID)%COUNT_FACETS
367              RADSQ = DES_RADIUS(LL)*DES_RADIUS(LL)
368     
369              particle_max(:) = des_pos_new(:, LL) + des_radius(LL)
370              particle_min(:) = des_pos_new(:, LL) - des_radius(LL)
371     
372              DO CELL_COUNT = 1, cellneighbor_facet_num(cell_id)
373     
374                 axis = cellneighbor_facet(cell_id)%extentdir(cell_count)
375     
376                 NF = cellneighbor_facet(cell_id)%p(cell_count)
377     
378     ! Compute particle-particle VDW cohesive short-range forces
379                 IF(USE_COHESION .AND. VAN_DER_WAALS) THEN
380     
381                    CALL ClosestPtPointTriangle(DES_POS_NEW(:,LL),          &
382                       VERTEX(:,:,NF), CLOSEST_PT(:))
383     
384                    DIST(:) = CLOSEST_PT(:) - DES_POS_NEW(:,LL)
385                    DISTSQ = DOT_PRODUCT(DIST, DIST)
386                    R_LM = 2*DES_RADIUS(LL)
387     
388                    IF(DISTSQ < (R_LM+WALL_VDW_OUTER_CUTOFF)**2) THEN
389                       IF(DISTSQ > (WALL_VDW_INNER_CUTOFF+R_LM)**2) THEN
390                          DistApart = (SQRT(DISTSQ)-R_LM)
391                          FORCE_COH = WALL_HAMAKER_CONSTANT * DES_RADIUS(LL) /   &
392                             (12d0*DistApart**2)*(Asperities/(Asperities +  &
393                             DES_RADIUS(LL)) + ONE/(ONE+Asperities/         &
394                             DistApart)**2)
395                       ELSE
396                          FORCE_COH = 2d0*PI*SURFACE_ENERGY*DES_RADIUS(LL)* &
397                             (Asperities/(Asperities+DES_RADIUS(LL)) + ONE/ &
398                             (ONE+Asperities/WALL_VDW_INNER_CUTOFF)**2 )
399                       ENDIF
400                       FC(:,LL) = FC(:,LL) + DIST(:)*FORCE_COH/SQRT(DISTSQ)
401                    ENDIF
402                 ENDIF
403     
404                 if (cellneighbor_facet(cell_id)%extentmin(cell_count) >    &
405                    particle_max(axis)) then
406                    call remove_collision(LL, nf, wall_collision_facet_id,wall_collision_PFT)
407                    cycle
408                 endif
409     
410                 if (cellneighbor_facet(cell_id)%extentmax(cell_count) <    &
411                    particle_min(axis)) then
412                    call remove_collision(LL, nf, wall_collision_facet_id,wall_collision_PFT)
413                    cycle
414                 endif
415     
416     ! Checking all the facets is time consuming due to the expensive
417     ! separating axis test. Remove this facet from contention based on
418     !a simple orthogonal projection test.
419     
420     ! Parametrize a line as p = p_0 + t normal and intersect with the
421     ! triangular plane. If t>0, then point is on the non-fluid side of
422     ! the plane. If the plane normal is assumed to point toward the fluid.
423     
424     ! -undefined, because non zero values will imply the sphere center
425     ! is on the non-fluid side of the plane. Since the testing
426     ! is with extended plane, this could very well happen even
427     ! when the particle is well inside the domain (assuming the plane
428     ! normal points toward the fluid). See the pic below. So check
429     ! only when line_t is negative
430     
431     !                 \   Solid  /
432     !                  \  Side  /
433     !                   \      /
434     !                    \    /
435     ! Wall 1, fluid side  \  /  Wall 2, fluid side
436     !                      \/
437     !                        o particle
438     !
439     ! line_t will be positive for wall 1 (incorrectly indicating center
440     ! is outside the domain) and line_t will be negative for wall 2.
441     !
442     ! Therefore, only stick with this test when line_t is negative and let
443     ! the separating axis test take care of the other cases.
444     
445     ! Since this is for checking static config, line's direction is the
446     ! same as plane's normal. For moving particles, the line's normal will
447     ! be along the point joining new and old positions.
448     
449                 line_t = DOT_PRODUCT(VERTEX(1,:,NF) - des_pos_new(:,LL),&
450                    NORM_FACE(:,NF))
451     
452     ! k - rad >= tol_orth, where k = -line_t, then orthogonal
453     ! projection is false. Substituting for k
454     ! => line_t + rad <= -tol_orth
455     ! choosing tol_orth = 0.01% of des_radius = 0.0001*des_radius
456     
457     ! Orthogonal projection will detect false positives even
458     ! when the particle does not overlap the triangle.
459     ! However, if the orthogonal projection shows no overlap, then
460     ! that is a big fat negative and overlaps are not possible.
461                 if((line_t.le.-1.0001d0*des_radius(LL))) then  ! no overlap
462                    call remove_collision(LL,nf,wall_collision_facet_id,wall_collision_PFT)
463                    CYCLE
464                 ENDIF
465     
466                 CALL ClosestPtPointTriangle(DES_POS_NEW(:,LL),             &
467                    VERTEX(:,:,NF), CLOSEST_PT(:))
468     
469                 DIST(:) = CLOSEST_PT(:) - DES_POS_NEW(:,LL)
470                 DISTSQ = DOT_PRODUCT(DIST, DIST)
471     
472                 IF(DISTSQ .GE. RADSQ) THEN !No overlap exists
473                    call remove_collision(LL,nf,wall_collision_facet_id,wall_collision_PFT)
474                    CYCLE
475                 ENDIF
476     
477                 MAX_DISTSQ = DISTSQ
478                 MAX_NF = NF
479     
480     ! Assign the collision normal based on the facet with the
481     ! largest overlap.
482                 NORMAL(:) = DIST(:)/sqrt(DISTSQ)
483     
484                       !NORMAL(:) = -NORM_FACE(:,MAX_NF)
485                    !facet's normal is correct normal only when the
486                    !intersection is with the face. When the intersection
487                    !is with edge or vertex, then the normal is
488                    !based on closest pt and sphere center. The
489                    !definition above of the normal is generic enough to
490                    !account for differences between vertex, edge, and facet.
491     
492     ! Calculate the particle/wall overlap.
493                    DISTMOD = SQRT(MAX_DISTSQ)
494                    OVERLAP_N = DES_RADIUS(LL) - DISTMOD
495     
496     ! Calculate the translational relative velocity for a contacting particle pair
497                    CALL CFRELVEL_WALL2(LL, V_REL_TRANS_NORM, &
498                       TANGENT, NORMAL, DISTMOD)
499     
500     ! Calculate the spring model parameters.
501                    phaseLL = PIJK(LL,5)
502                    ! Hertz vs linear spring-dashpot contact model
503                    IF (DES_COLL_MODEL_ENUM .EQ. HERTZIAN) THEN
504                       sqrt_overlap = SQRT(OVERLAP_N)
505                       KN_DES_W = hert_kwn(phaseLL)*sqrt_overlap
506                       KT_DES_W = hert_kwt(phaseLL)*sqrt_overlap
507                       sqrt_overlap = SQRT(sqrt_overlap)
508                       ETAN_DES_W = DES_ETAN_WALL(phaseLL)*sqrt_overlap
509                       ETAT_DES_W = DES_ETAT_WALL(phaseLL)*sqrt_overlap
510                    ELSE
511                       KN_DES_W = KN_W
512                       KT_DES_W = KT_W
513                       ETAN_DES_W = DES_ETAN_WALL(phaseLL)
514                       ETAT_DES_W = DES_ETAT_WALL(phaseLL)
515                    ENDIF
516     
517     ! Calculate the normal contact force
518                    FNS1(:) = -KN_DES_W * OVERLAP_N * NORMAL(:)
519                    FNS2(:) = -ETAN_DES_W * V_REL_TRANS_NORM * NORMAL(:)
520                    FNORM(:) = FNS1(:) + FNS2(:)
521     
522     ! Calculate the tangential displacement. Note that only the maximum
523     ! wall collision is considered. Therefore, enduring contact can exist
524     ! from facet-to-facet as long as the particle remains in contact with
525     ! one or more facets.
526     
527                    pft = get_collision(LL,nf,wall_collision_facet_id,wall_collision_PFT)
528     
529                    IF(sum(abs(PFT(:))) .gt. small_number) THEN
530                       OVERLAP_T(:) = TANGENT(:) * DTSOLID
531                    ELSE
532                       IF(V_REL_TRANS_NORM .GT. ZERO) THEN
533                          DTSOLID_TMP = OVERLAP_N/(V_REL_TRANS_NORM)
534                       ELSEIF(V_REL_TRANS_NORM .LT. ZERO) THEN
535                          DTSOLID_TMP = DTSOLID
536                       ELSE
537                          DTSOLID_TMP = OVERLAP_N /                         &
538                             (V_REL_TRANS_NORM+SMALL_NUMBER)
539                       ENDIF
540                       OVERLAP_T(:) = TANGENT(:) * MIN(DTSOLID,DTSOLID_TMP)
541                    ENDIF
542     
543     ! Update the tangential history.
544                    PFT(:) = PFT(:) + OVERLAP_T(:)
545     
546                    FORCE_HISTORY(:) = PFT(:) - &
547                       DOT_PRODUCT(PFT(:),NORMAL)*NORMAL(:)
548     
549     ! Calculate the tangential collision force.
550                    FTS1(:) = -KT_DES_W * FORCE_HISTORY(:)
551                    FTS2(:) = -ETAT_DES_W * TANGENT(:)
552                    FTAN(:) =  FTS1(:) + FTS2(:)
553     
554                    PARTICLE_SLIDE = .FALSE.
555     
556     ! Check for Coulombs friction law and limit the maximum value of the
557     ! tangential force on a particle in contact with a wall.
558                    FTMD = DOT_PRODUCT(FTAN, FTAN)
559                    FNMD = DOT_PRODUCT(FNORM,FNORM)
560                    IF (FTMD.GT.(MEW_W*MEW_W*FNMD)) THEN
561                       PARTICLE_SLIDE = .TRUE.
562                       IF(all(TANGENT.EQ.zero)) THEN
563                          FTAN(:) =  MEW_W * FTAN(:) * SQRT(FNMD/FTMD)
564                       ELSE
565                          FTAN(:) = -MEW_W * TANGENT(:) * SQRT(FNMD/dot_product(TANGENT,TANGENT))
566                       ENDIF
567                    ENDIF
568     
569     ! Add the collision force to the total forces acting on the particle.
570                    FC(:,LL) = FC(:,LL) + FNORM(:) + FTAN(:)
571     
572     ! Add the torque: The particle radius is used as the moment arm
573                    CROSSP = DES_CROSSPRDCT(NORMAL, FTAN)
574                    TOW(:,LL) = TOW(:,LL) + DES_RADIUS(LL)*CROSSP(:)
575     
576     ! Save the tangential displacement history with the correction of Coulomb's law
577                    IF (PARTICLE_SLIDE) THEN
578     ! Since FT might be corrected during the call to cfslide, the tangential
579     ! displacement history needs to be changed accordingly
580                       PFT(:) = -( FTAN(:) - FTS2(:) ) / KT_DES_W
581                    ELSE
582                       PFT(:) = FORCE_HISTORY(:)
583                    ENDIF
584     
585                    call update_collision(pft,LL,nf,wall_collision_facet_id,wall_collision_PFT)
586     
587              ENDDO
588     
589     !         ENDIF
590           ENDDO
591     !$omp end do
592     !$omp end parallel
593     
594           RETURN
595     
596            contains
597     
598              function get_collision(LLL,facet_id,wall_collision_facet_id,wall_collision_PFT)
599                use error_manager
600                implicit none
601                double precision, dimension(DIMN) :: get_collision
602                Integer, intent(in) :: LLL,facet_id
603                INTEGER, DIMENSION(:,:), intent(inout) :: wall_collision_facet_id
604                DOUBLE PRECISION, DIMENSION(:,:,:), intent(inout) :: wall_collision_PFT
605                integer :: cc, free_index
606     
607                free_index = -1
608     
609                do cc = 1, COLLISION_ARRAY_MAX
610                   if (facet_id == wall_collision_facet_id(cc,LLL)) then
611                      get_collision(:) = wall_collision_PFT(:,cc,LLL)
612                      return
613                   else if (-1 == wall_collision_facet_id(cc,LLL)) then
614                      free_index = cc
615                   endif
616                enddo
617                if(-1 == free_index) then
618                   CALL INIT_ERR_MSG("CALC_COLLISION_WALL_MOD: GET_COLLISION")
619                   WRITE(ERR_MSG, 1100)
620                   CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
621                else
622                   wall_collision_facet_id(free_index,LLL) = facet_id
623                   wall_collision_PFT(:,free_index,LLL) = ZERO
624                   get_collision(:) = wall_collision_PFT(:,free_index,LLL)
625                   return
626                endif
627     
628     1100       FORMAT('Error: COLLISION_ARRAY_MAX too small. ')
629     
630              end function get_collision
631     
632              subroutine update_collision(pft,LLL,facet_id,wall_collision_facet_id,wall_collision_PFT)
633                use error_manager
634                implicit none
635                double precision, dimension(DIMN), intent(in) :: pft
636                Integer, intent(in) :: LLL,facet_id
637                INTEGER, DIMENSION(:,:), intent(inout) :: wall_collision_facet_id
638                DOUBLE PRECISION, DIMENSION(:,:,:), intent(inout) :: wall_collision_PFT
639                integer :: cc, free_index
640     
641                do cc = 1, COLLISION_ARRAY_MAX
642                   if (facet_id == wall_collision_facet_id(cc,LLL)) then
643                      wall_collision_PFT(:,cc,LLL) = PFT(:)
644                      return
645                   endif
646                enddo
647     
648                CALL INIT_ERR_MSG("CALC_COLLISION_WALL_MOD: UPDATE_COLLISION")
649                WRITE(ERR_MSG, 1100)
650                CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
651     
652     1100       FORMAT('Error: COLLISION_ARRAY_MAX too small. ')
653     
654              end subroutine update_collision
655     
656              subroutine remove_collision(LLL,facet_id,wall_collision_facet_id,wall_collision_PFT)
657                use error_manager
658                implicit none
659                Integer, intent(in) :: LLL,facet_id
660                INTEGER, DIMENSION(:,:), intent(inout) :: wall_collision_facet_id
661                DOUBLE PRECISION, DIMENSION(:,:,:), intent(inout) :: wall_collision_PFT
662                integer :: cc
663     
664                do cc = 1, COLLISION_ARRAY_MAX
665                   if (facet_id == wall_collision_facet_id(cc,LLL)) then
666                      wall_collision_facet_id(cc,LLL) = -1
667                      return
668                   endif
669                enddo
670     
671     1100       FORMAT('Error: COLLISION_ARRAY_MAX too small. ')
672     
673              end subroutine remove_collision
674     
675         END SUBROUTINE CALC_DEM_FORCE_WITH_WALL_STL
676     
677     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
678     !                                                                      !
679     !                                                                      !
680     !                                                                      !
681     !                                                                      !
682     !                                                                      !
683     !                                                                      !
684     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
685         SUBROUTINE write_this_facet_and_part(FID, PID)
686           USE run
687           USE param1
688           USE discretelement
689           USE geometry
690           USE compar
691           USE constant
692           USE cutcell
693           USE funits
694           USE indices
695           USE physprop
696           USE parallel
697           USE stl
698           USE des_stl_functions
699           Implicit none
700           !facet id and particle id
701           Integer, intent(in) :: fid, pid
702           Integer :: stl_unit, vtp_unit , k
703           CHARACTER(LEN=100) :: stl_fname, vtp_fname
704           real :: temp_array(3)
705     
706           stl_unit = 1001
707           vtp_unit = 1002
708     
709           WRITE(vtp_fname,'(A,"_OFFENDING_PARTICLE",".vtp")') TRIM(RUN_NAME)
710           WRITE(stl_fname,'(A,"_STL_FACE",".stl")') TRIM(RUN_NAME)
711     
712           open(vtp_unit, file = vtp_fname, form='formatted')
713           open(stl_unit, file = stl_fname, form='formatted')
714     
715           write(vtp_unit,"(a)") '<?xml version="1.0"?>'
716           write(vtp_unit,"(a,es24.16,a)") '<!-- time =',s_time,'s -->'
717           write(vtp_unit,"(a,a)") '<VTKFile type="PolyData"',&
718                ' version="0.1" byte_order="LittleEndian">'
719           write(vtp_unit,"(3x,a)") '<PolyData>'
720           write(vtp_unit,"(6x,a,i10.10,a,a)")&
721                '<Piece NumberOfPoints="',1,'" NumberOfVerts="0" ',&
722                'NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="0">'
723           write(vtp_unit,"(9x,a)")&
724                '<PointData Scalars="Diameter" Vectors="Velocity">'
725           write(vtp_unit,"(12x,a)")&
726                '<DataArray type="Float32" Name="Diameter" format="ascii">'
727           write (vtp_unit,"(15x,es13.6)") (2*des_radius(pid))
728           write(vtp_unit,"(12x,a)") '</DataArray>'
729     
730           temp_array = zero
731           temp_array(:) = des_vel_new(:,pid)
732           write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
733                'Name="Velocity" NumberOfComponents="3" format="ascii">'
734           write (vtp_unit,"(15x,3(es13.6,3x))")&
735                ((temp_array(k)),k=1,3)
736           write(vtp_unit,"(12x,a,/9x,a)") '</DataArray>','</PointData>'
737           ! skip cell data
738           write(vtp_unit,"(9x,a)") '<CellData></CellData>'
739     
740           temp_array = zero
741           temp_array(1:dimn) = des_pos_new(1:dimn, pid)
742           write(vtp_unit,"(9x,a)") '<Points>'
743           write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
744                'Name="Position" NumberOfComponents="3" format="ascii">'
745           write (vtp_unit,"(15x,3(es13.6,3x))")&
746                ((temp_array(k)),k=1,3)
747           write(vtp_unit,"(12x,a,/9x,a)")'</DataArray>','</Points>'
748           ! Write tags for data not included (vtp format style)
749           write(vtp_unit,"(9x,a,/9x,a,/9x,a,/9x,a)")'<Verts></Verts>',&
750                '<Lines></Lines>','<Strips></Strips>','<Polys></Polys>'
751           write(vtp_unit,"(6x,a,/3x,a,/a)")&
752                '</Piece>','</PolyData>','</VTKFile>'
753     
754           !Now write the facet info
755     
756           write(stl_unit,*)'solid vcg'
757     
758           write(stl_unit,*) '   facet normal ', NORM_FACE(1:3,FID)
759           write(stl_unit,*) '      outer loop'
760           write(stl_unit,*) '         vertex ', VERTEX(1,1:3,FID)
761           write(stl_unit,*) '         vertex ', VERTEX(2,1:3,FID)
762           write(stl_unit,*) '         vertex ', VERTEX(3,1:3,FID)
763           write(stl_unit,*) '      endloop'
764           write(stl_unit,*) '   endfacet'
765     
766           write(stl_unit,*)'endsolid vcg'
767     
768           close(vtp_unit, status = 'keep')
769           close(stl_unit, status = 'keep')
770     
771         end SUBROUTINE write_this_facet_and_part
772     
773     !----------------------------------------------------------------------!
774     !                                                                      !
775     !                                                                      !
776     !                                                                      !
777     !                                                                      !
778     !                                                                      !
779     !----------------------------------------------------------------------!
780           SUBROUTINE CFRELVEL_WALL2(LL,  VRN, VSLIP, NORM, DIST_LI)
781     
782           USE discretelement
783           USE param1
784     
785           IMPLICIT NONE
786     
787           INTEGER, INTENT(IN) :: LL
788           DOUBLE PRECISION, DIMENSION(DIMN), INTENT(IN) :: NORM
789           DOUBLE PRECISION, DIMENSION(DIMN), INTENT(OUT):: VSLIP
790           DOUBLE PRECISION, INTENT(OUT):: VRN
791     ! distance between particles
792           DOUBLE PRECISION, INTENT(IN) :: DIST_LI
793     
794     !-----------------------------------------------
795     ! Local variables
796     !-----------------------------------------------
797     
798           DOUBLE PRECISION, DIMENSION(DIMN) :: V_ROT, OMEGA_SUM, VRELTRANS
799     
800     ! distance from the contact point to the particle centers
801           DOUBLE PRECISION DIST_CL
802     
803     ! translational relative velocity
804           VRELTRANS(:) = DES_VEL_NEW(:,LL)
805     
806     ! rotational contribution  : v_rot
807     ! calculate the distance from the particle center to the wall
808           DIST_CL = DIST_LI         !- DES_RADIUS(L)
809           OMEGA_SUM(:) = OMEGA_NEW(:,LL)*DIST_CL
810     
811           V_ROT = DES_CROSSPRDCT(OMEGA_SUM, NORM)
812     
813     ! total relative velocity
814           VRELTRANS(:) =  DES_VEL_NEW(:,LL) + V_ROT(:)
815     
816     ! normal component of relative velocity (scalar)
817           VRN = DOT_PRODUCT(VRELTRANS,NORM)
818     
819     ! slip velocity of the contact point
820     ! Equation (8) in Tsuji et al. 1992
821           VSLIP(:) =  VRELTRANS(:) - VRN*NORM(:)
822     
823           RETURN
824           END SUBROUTINE CFRELVEL_WALL2
825     
826         end module CALC_COLLISION_WALL
827     
828