File: RELATIVE:/../../../mfix.git/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     
18           CONTAINS
19     
20     !VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV!
21     !                                                                      !
22     !                                                                      !
23     !                                                                      !
24     !                                                                      !
25     !                                                                      !
26     !                                                                      !
27     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
28           SUBROUTINE CALC_DEM_FORCE_WITH_WALL_STL
29     
30           USE run
31           USE param1
32           USE desgrid
33           USE discretelement
34           USE geometry
35           USE compar
36           USE constant
37           USE indices
38           USE stl
39           USE des_stl_functions
40           USE functions
41           Implicit none
42     
43           INTEGER :: LL
44           INTEGER IJK, NF
45           DOUBLE PRECISION OVERLAP_N, SQRT_OVERLAP
46     
47           DOUBLE PRECISION V_REL_TRANS_NORM, DISTSQ, RADSQ, CLOSEST_PT(DIMN)
48     ! local normal and tangential forces
49           DOUBLE PRECISION NORMAL(DIMN), VREL_T(DIMN), DIST(DIMN), DISTMOD
50           DOUBLE PRECISION, DIMENSION(DIMN) :: FTAN, FNORM, OVERLAP_T
51     
52           LOGICAL :: DES_LOC_DEBUG
53           INTEGER :: CELL_ID, cell_count
54           INTEGER :: PHASELL
55     
56           DOUBLE PRECISION :: TANGENT(DIMN)
57           DOUBLE PRECISION :: FTMD, FNMD
58     ! local values used spring constants and damping coefficients
59           DOUBLE PRECISION ETAN_DES_W, ETAT_DES_W, KN_DES_W, KT_DES_W
60     
61           double precision :: MAG_OVERLAP_T
62     
63           double precision :: line_t
64     ! flag to tell if the orthogonal projection of sphere center to
65     ! extended plane detects an overlap
66     
67           DOUBLE PRECISION :: MAX_DISTSQ, DISTAPART, FORCE_COH, R_LM
68           INTEGER :: MAX_NF, axis
69           DOUBLE PRECISION, DIMENSION(3) :: PARTICLE_MIN, PARTICLE_MAX
70     
71           DES_LOC_DEBUG = .false. ;      DEBUG_DES = .false.
72           FOCUS_PARTICLE = -1
73     
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)         &
81     !$omp shared(max_pip,focus_particle,debug_des,no_neighboring_facet_des,&
82     !$omp    pijk,dg_pijk,list_facet_at_des,i_of,j_of,k_of,des_pos_new,    &
83     !$omp    des_radius,cellneighbor_facet_num,cellneighbor_facet,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 (NO_NEIGHBORING_FACET_DES(DG_PIJK(LL))) THEN
103              IF(CELLNEIGHBOR_FACET_NUM(CELL_ID) < 1) THEN
104                 WALL_COLLISION_FACET_ID(:,LL) = -1
105                 WALL_COLLISION_PFT(:,:,LL) = 0.0d0
106                 CYCLE
107              ENDIF
108     
109     ! Check particle LL for wall contacts
110              RADSQ = DES_RADIUS(LL)*DES_RADIUS(LL)
111     
112              particle_max(:) = des_pos_new(:, LL) + des_radius(LL)
113              particle_min(:) = des_pos_new(:, LL) - des_radius(LL)
114     
115              DO CELL_COUNT = 1, cellneighbor_facet_num(cell_id)
116     
117                 axis = cellneighbor_facet(cell_id)%extentdir(cell_count)
118     
119                 NF = cellneighbor_facet(cell_id)%p(cell_count)
120     
121     ! Compute particle-particle VDW cohesive short-range forces
122                 IF(USE_COHESION .AND. VAN_DER_WAALS) THEN
123     
124                    CALL ClosestPtPointTriangle(DES_POS_NEW(:,LL),          &
125                       VERTEX(:,:,NF), CLOSEST_PT(:))
126     
127                    DIST(:) = CLOSEST_PT(:) - DES_POS_NEW(:,LL)
128                    DISTSQ = DOT_PRODUCT(DIST, DIST)
129                    R_LM = 2*DES_RADIUS(LL)
130     
131                    IF(DISTSQ < (R_LM+WALL_VDW_OUTER_CUTOFF)**2) THEN
132                       IF(DISTSQ > (WALL_VDW_INNER_CUTOFF+R_LM)**2) THEN
133                          DistApart = (SQRT(DISTSQ)-R_LM)
134                          FORCE_COH = WALL_HAMAKER_CONSTANT*DES_RADIUS(LL) /&
135                             (12d0*DistApart**2)*(Asperities/(Asperities +  &
136                             DES_RADIUS(LL)) + ONE/(ONE+Asperities/         &
137                             DistApart)**2)
138                       ELSE
139                          FORCE_COH = 2d0*PI*SURFACE_ENERGY*DES_RADIUS(LL)* &
140                             (Asperities/(Asperities+DES_RADIUS(LL)) + ONE/ &
141                             (ONE+Asperities/WALL_VDW_INNER_CUTOFF)**2 )
142                       ENDIF
143                       FC(:,LL) = FC(:,LL) + DIST(:)*FORCE_COH/SQRT(DISTSQ)
144                    ENDIF
145                 ENDIF
146     
147                 if (cellneighbor_facet(cell_id)%extentmin(cell_count) >    &
148                    particle_max(axis)) then
149                    call remove_collision(LL, nf, wall_collision_facet_id)
150                    cycle
151                 endif
152     
153                 if (cellneighbor_facet(cell_id)%extentmax(cell_count) <    &
154                    particle_min(axis)) then
155                    call remove_collision(LL, nf, wall_collision_facet_id)
156                    cycle
157                 endif
158     
159     ! Checking all the facets is time consuming due to the expensive
160     ! separating axis test. Remove this facet from contention based on
161     ! a simple orthogonal projection test.
162     
163     ! Parametrize a line as p = p_0 + t normal and intersect with the
164     ! triangular plane. If t>0, then point is on the non-fluid side of
165     ! the plane. If the plane normal is assumed to point toward the fluid.
166     
167     ! -undefined, because non zero values will imply the sphere center
168     ! is on the non-fluid side of the plane. Since the testing
169     ! is with extended plane, this could very well happen even
170     ! when the particle is well inside the domain (assuming the plane
171     ! normal points toward the fluid). See the pic below. So check
172     ! only when line_t is negative
173     
174     !                 \   Solid  /
175     !                  \  Side  /
176     !                   \      /
177     !                    \    /
178     ! Wall 1, fluid side  \  /  Wall 2, fluid side
179     !                      \/
180     !                        o particle
181     !
182     ! line_t will be positive for wall 1 (incorrectly indicating center
183     ! is outside the domain) and line_t will be negative for wall 2.
184     !
185     ! Therefore, only stick with this test when line_t is negative and let
186     ! the separating axis test take care of the other cases.
187     
188     ! Since this is for checking static config, line's direction is the
189     ! same as plane's normal. For moving particles, the line's normal will
190     ! be along the point joining new and old positions.
191     
192                 line_t = DOT_PRODUCT(VERTEX(1,:,NF) - des_pos_new(:,LL),&
193                    NORM_FACE(:,NF))
194     
195     ! k - rad >= tol_orth, where k = -line_t, then orthogonal
196     ! projection is false. Substituting for k
197     ! => line_t + rad <= -tol_orth
198     ! choosing tol_orth = 0.01% of des_radius = 0.0001*des_radius
199     
200     ! Orthogonal projection will detect false positives even
201     ! when the particle does not overlap the triangle.
202     ! However, if the orthogonal projection shows no overlap, then
203     ! that is a big fat negative and overlaps are not possible.
204                 if((line_t.le.-1.0001d0*des_radius(LL))) then  ! no overlap
205                    call remove_collision(LL,nf,wall_collision_facet_id)
206                    CYCLE
207                 ENDIF
208     
209                 CALL ClosestPtPointTriangle(DES_POS_NEW(:,LL),             &
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) 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           use stl_preproc_des
317           use error_manager
318     
319           IMPLICIT NONE
320     
321           DOUBLE PRECISION :: GET_COLLISION(DIMN)
322           INTEGER, INTENT(IN) :: LLL,FACET_ID
323           INTEGER, INTENT(INOUT) :: WALL_COLLISION_FACET_ID(:,:)
324           DOUBLE PRECISION, INTENT(INOUT) :: WALL_COLLISION_PFT(:,:,:)
325           INTEGER :: CC, FREE_INDEX, LC, dgIJK
326     
327           free_index = -1
328     
329           do cc = 1, COLLISION_ARRAY_MAX
330              if (facet_id == wall_collision_facet_id(cc,LLL)) then
331                 get_collision(:) = wall_collision_PFT(:,cc,LLL)
332                 return
333              else if (-1 == wall_collision_facet_id(cc,LLL)) then
334                 free_index = cc
335              endif
336           enddo
337     
338     ! Overwrite old data. This is needed because a particle moving from
339     ! one dg cell to another may no longer 'see' an STL before it moved
340     ! out of contact range. Therefore, the 'remove_collision' function
341     ! does not get called to cleanup the stale data.
342           if(-1 == free_index) then
343              dgIJK=DG_PIJK(LLL)
344              cc_lp: do cc=1, COLLISION_ARRAY_MAX
345                 do lc=1, cellneighbor_facet_num(dgIJK)
346                    if(wall_collision_facet_id(cc,LLL) == &
347                       cellneighbor_facet(dgIJK)%p(LC))  cycle cc_lp
348                 enddo
349                 free_index = cc
350                 exit cc_lp
351              enddo cc_lp
352           endif
353     
354     
355           if(-1 == free_index) then
356     
357              do cc = 1, COLLISION_ARRAY_MAX
358                 call write_this_stl(wall_collision_facet_id(cc,LLL))
359              enddo
360              call write_stls_this_dg(dg_pijk(LLL))
361     
362              CALL INIT_ERR_MSG("CALC_COLLISION_WALL_MOD: GET_COLLISION")
363              WRITE(ERR_MSG, 1100) LLL, CC
364              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
365     
366           else
367              wall_collision_facet_id(free_index,LLL) = facet_id
368              wall_collision_PFT(:,free_index,LLL) = ZERO
369              get_collision(:) = wall_collision_PFT(:,free_index,LLL)
370              return
371           endif
372     
373      1100 FORMAT('Error: COLLISION_ARRAY_MAX too small.'/'Particle: ',I9,/ &
374                'Facet:    ',I9)
375     
376           END FUNCTION GET_COLLISION
377     
378     
379     !......................................................................!
380     !  Function: UPDATE_COLLISION                                          !
381     !                                                                      !
382     !  Purpose: Update the integrated (t0->t) tangential displacement.     !
383     !......................................................................!
384           SUBROUTINE UPDATE_COLLISION(PFT, LLL, FACET_ID,                  &
385              WALL_COLLISION_FACET_ID, WALL_COLLISION_PFT)
386     
387           use error_manager
388           implicit none
389     
390           DOUBLE PRECISION, INTENT(IN) :: PFT(DIMN)
391           INTEGER, INTENT(IN) :: LLL,FACET_ID
392           INTEGER, INTENT(IN) :: WALL_COLLISION_FACET_ID(:,:)
393           DOUBLE PRECISION, INTENT(INOUT) :: WALL_COLLISION_PFT(:,:,:)
394           INTEGER :: CC
395     
396           do cc = 1, COLLISION_ARRAY_MAX
397              if (facet_id == wall_collision_facet_id(cc,LLL)) then
398                 wall_collision_PFT(:,cc,LLL) = PFT(:)
399                 return
400              endif
401           enddo
402     
403           CALL INIT_ERR_MSG("CALC_COLLISION_WALL_MOD: UPDATE_COLLISION")
404           WRITE(ERR_MSG, 1100)
405           CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
406     
407      1100 FORMAT('Error: COLLISION_ARRAY_MAX too small. ')
408     
409           END SUBROUTINE UPDATE_COLLISION
410     
411     !......................................................................!
412     !  Function: REMOVE_COLLISION                                          !
413     !                                                                      !
414     !  Purpose: Clear the integrated (t0->t) tangential displacement once  !
415     !  the collision is over (contact ended).                              !
416     !......................................................................!
417           SUBROUTINE REMOVE_COLLISION(LLL,FACET_ID,WALL_COLLISION_FACET_ID)
418     
419           use error_manager
420     
421           IMPLICIT NONE
422     
423           INTEGER, INTENT(IN) :: LLL,FACET_ID
424           INTEGER, INTENT(INOUT) :: WALL_COLLISION_FACET_ID(:,:)
425           INTEGER :: CC
426     
427           DO CC = 1, COLLISION_ARRAY_MAX
428              IF (FACET_ID == WALL_COLLISION_FACET_ID(CC,LLL)) THEN
429                 WALL_COLLISION_FACET_ID(CC,LLL) = -1
430                 RETURN
431              ENDIF
432           ENDDO
433     
434           END SUBROUTINE REMOVE_COLLISION
435     
436           END SUBROUTINE CALC_DEM_FORCE_WITH_WALL_STL
437     
438     
439     !----------------------------------------------------------------------!
440     !                                                                      !
441     !  Subroutine: CFRELVEL_WALL                                           !
442     !                                                                      !
443     !  Purpose: Calculate the normal and tangential components of the      !
444     !  relative velocity between a particle and wall contact.              !
445     !                                                                      !
446     !  Comments: Only the magnitude of the normal component is returned    !
447     !  whereas the full tangential vector is returned.                     !
448     !----------------------------------------------------------------------!
449           SUBROUTINE CFRELVEL_WALL(LL, VRN, VRT, NORM, DIST)
450     
451     ! Particle translational velocity
452           use discretelement, only: DES_VEL_NEW
453     ! Particle rotational velocity
454           use discretelement, only: OMEGA_NEW
455     ! Spatial array size (parameter)
456           use discretelement, only: DIMN
457     ! Function for calculating the cross prodcut
458           use discretelement, only: DES_CROSSPRDCT
459     
460           IMPLICIT NONE
461     
462     ! Dummy arguments:
463     !---------------------------------------------------------------------//
464     ! Particle index.
465           INTEGER, INTENT(IN) :: LL
466     ! Magnitude of the total relative translational velocity.
467           DOUBLE PRECISION, INTENT(OUT):: VRN
468     ! Total relative translational velocity (vector).
469           DOUBLE PRECISION, DIMENSION(DIMN), INTENT(OUT):: VRT
470     ! Unit normal from particle center to closest point on stl (wall)
471           DOUBLE PRECISION, DIMENSION(DIMN), INTENT(IN) :: NORM
472     ! Distance between particle center and stl (wall).
473           DOUBLE PRECISION, INTENT(IN) :: DIST
474     
475     ! Local variables
476     !---------------------------------------------------------------------//
477     ! Additional relative translational motion due to rotation
478           DOUBLE PRECISION, DIMENSION(DIMN) :: V_ROT
479     ! Total relative velocity at contact point
480           DOUBLE PRECISION, DIMENSION(DIMN) :: VRELTRANS
481     
482     ! Total relative velocity + rotational contribution
483           V_ROT = DIST*OMEGA_NEW(:,LL)
484           VRELTRANS(:) =  DES_VEL_NEW(:,LL) + DES_CROSSPRDCT(V_ROT, NORM)
485     
486     ! magnitude of normal component of relative velocity (scalar)
487           VRN = DOT_PRODUCT(VRELTRANS,NORM)
488     
489     ! total relative translational slip velocity at the contact point
490     ! Equation (8) in Tsuji et al. 1992
491           VRT(:) =  VRELTRANS(:) - VRN*NORM(:)
492     
493           RETURN
494           END SUBROUTINE CFRELVEL_WALL
495     
496           end module CALC_COLLISION_WALL
497     
498