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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: CALC_FORCE_DEM                                          !
4     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  !
5     !                                                                      !
6     !  Purpose: Calculate contact force and torque on particle from        !
7     !           particle-particle and particle-wall collisions. Treats     !
8     !           wall interaction also as a two-particle interaction but    !
9     !           accounting for the wall properties                         !
10     !                                                                      !
11     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
12     #include "version.inc"
13           SUBROUTINE CALC_FORCE_DEM
14     
15     ! Modules
16     !---------------------------------------------------------------------//
17           USE calc_collision_wall
18           USE constant, ONLY: Pi
19           USE derived_types, only: multisap, boxhandle
20           USE des_thermo
21           USE des_thermo_cond
22           USE discretelement
23           USE run
24           use pair_manager
25           use param1, only: one, small_number, zero
26     
27           IMPLICIT NONE
28     
29     ! Local variables
30     !---------------------------------------------------------------------//
31     ! percent of particle radius when excess overlap will be flagged
32           DOUBLE PRECISION, PARAMETER :: flag_overlap = 0.20d0
33     ! particle no. indices
34           INTEGER :: I, LL, cc, CC_START, CC_END
35     ! the overlap occuring between particle-particle or particle-wall
36     ! collision in the normal direction
37           DOUBLE PRECISION :: OVERLAP_N, OVERLAP_T(3)
38     ! square root of the overlap
39           DOUBLE PRECISION :: SQRT_OVERLAP
40     ! distance vector between two particle centers or between a particle
41     ! center and wall when the two surfaces are just at contact (i.e. no
42     ! overlap)
43           DOUBLE PRECISION :: R_LM,DIST_CI,DIST_CL
44     ! the normal and tangential components of the translational relative
45     ! velocity
46           DOUBLE PRECISION :: V_REL_TRANS_NORM, rad
47     ! distance vector between two particle centers or between a particle
48     ! center and wall at current and previous time steps
49           DOUBLE PRECISION :: DIST(3), NORMAL(3), DIST_MAG, POS(3)
50     ! tangent to the plane of contact at current time step
51           DOUBLE PRECISION :: VREL_T(3)
52     ! normal and tangential forces
53           DOUBLE PRECISION :: FN(3), FT(3)
54     ! temporary storage of force
55           DOUBLE PRECISION :: FC_TMP(3)
56     ! temporary storage of force for torque
57           DOUBLE PRECISION :: TOW_FORCE(3)
58     ! temporary storage of torque
59           DOUBLE PRECISION :: TOW_TMP(3,2)
60     ! temporary storage of conduction/radiation
61           DOUBLE PRECISION :: QQ_TMP
62     
63     ! store solids phase index of particle (i.e. pijk(np,5))
64           INTEGER :: PHASEI, PHASELL
65     ! local values used spring constants and damping coefficients
66           DOUBLE PRECISION :: ETAN_DES, ETAT_DES
67           DOUBLE PRECISION :: KN_DES, KT_DES
68     ! local values used for calculating cohesive forces
69           DOUBLE PRECISION :: FORCE_COH, EQ_RADIUS, DistApart
70     
71           LOGICAL, PARAMETER :: report_excess_overlap = .FALSE.
72     
73           DOUBLE PRECISION :: FNMD, FTMD, MAG_OVERLAP_T, TANGENT(3)
74     
75           integer :: nn, mm, box_id, box_id2
76           logical :: found
77           integer :: pair(2)
78     
79     !......................................................................!
80     
81     ! Initialize cohesive forces
82           IF(USE_COHESION) PostCohesive(:) = ZERO
83     
84           CALL CALC_DEM_FORCE_WITH_WALL_STL
85     
86     #ifdef do_sap
87           ! do nn=0, size(multisap%saps)-1
88           !    !print *,"nn = ",nn
89           !    if (.not.check_boxes(multisap%saps(nn))) ERROR_STOP __LINE__
90           !    if (.not.check_sort(multisap%saps(nn))) ERROR_STOP __LINE__
91           ! enddo
92     
93     !print *,"CALC_FORCE_DEM =================================================================================="
94     
95     print *," TOTAL NUM OF NEIGHBORS IS ",NEIGHBOR_INDEX(MAX_PIP-1)
96     
97     open (unit=123,file="neighbors.txt",action="write",status="replace")
98     
99     DO LL = 1, MAX_PIP
100        CC_START = 1
101        IF (LL.gt.1) CC_START = NEIGHBOR_INDEX(LL-1)
102        CC_END   = NEIGHBOR_INDEX(LL)
103     
104        DO CC = CC_START, CC_END-1
105           I  = NEIGHBORS(CC)
106           write (123,*) ll,i
107        enddo
108     enddo
109     
110     close (unit=123)
111     
112                  call reset_pairs(multisap%hashtable)
113                  do
114                     call get_pair(multisap%hashtable,pair)
115                     if (pair(1).eq.0 .and. pair(2).eq.0) exit
116     
117                     if ( des_radius(pair(1))+des_radius(pair(2))> sqrt(dot_product(DES_POS_NEW(pair(1),:)-DES_POS_NEW(pair(2),:),DES_POS_NEW(pair(1),:)-DES_POS_NEW(pair(2),:)))) then
118                        print *,"invalid pair: ",pair(1),pair(2)
119                        stop __LINE__
120                     endif
121                  enddo
122     #endif
123     
124     ! Check particle LL neighbor contacts
125     !---------------------------------------------------------------------//
126     
127     !!$omp parallel default(none) private(pos,rad,cc,cc_start,cc_end,ll,i,  &
128     !!$omp    overlap_n,vrel_t,v_rel_trans_norm,sqrt_overlap,dist,r_lm,     &
129     !!$omp    kn_des,kt_des,hert_kn,hert_kt,phasell,phasei,etan_des,        &
130     !!$omp    etat_des,fn,ft,overlap_t,tangent,mag_overlap_t,               &
131     !!$omp    eq_radius,distapart,force_coh,dist_mag,NORMAL,ftmd,fnmd,      &
132     !!$omp    dist_cl, dist_ci, fc_tmp, tow_tmp, tow_force, qq_tmp, box_id, box_id2, found)         &
133     !!$omp shared(max_pip,neighbors,neighbor_index,des_pos_new,des_radius,  &
134     !!$omp    des_coll_model_enum,kn,kt,pft_neighbor,pijk,                  &
135     !!$omp    des_etan,des_etat,mew,use_cohesion, calc_cond_des, dtsolid,   &
136     !!$omp    van_der_waals,vdw_outer_cutoff,vdw_inner_cutoff,              &
137     !!$omp    hamaker_constant,asperities,surface_energy,                   &
138     !!$omp    tow, fc, energy_eq, grav_mag, postcohesive, pmass, q_source, multisap, boxhandle)
139     
140     !!$omp do
141     
142           DO LL = 1, MAX_PIP
143              IF(IS_NONEXISTENT(LL)) CYCLE
144              pos = DES_POS_NEW(LL,:)
145              rad = DES_RADIUS(LL)
146     
147              CC_START = 1
148              IF (LL.gt.1) CC_START = NEIGHBOR_INDEX(LL-1)
149              CC_END   = NEIGHBOR_INDEX(LL)
150     
151              DO CC = CC_START, CC_END-1
152                 I  = NEIGHBORS(CC)
153                 IF(IS_NONEXISTENT(I)) CYCLE
154     
155                 R_LM = rad + DES_RADIUS(I)
156                 DIST(:) = DES_POS_NEW(I,:) - POS(:)
157                 DIST_MAG = dot_product(DIST,DIST)
158     
159                 FC_TMP(:) = ZERO
160     
161     ! Compute particle-particle VDW cohesive short-range forces
162                 IF(USE_COHESION .AND. VAN_DER_WAALS) THEN
163                    IF(DIST_MAG < (R_LM+VDW_OUTER_CUTOFF)**2) THEN
164                       EQ_RADIUS = 2d0 * DES_RADIUS(LL)*DES_RADIUS(I) /     &
165                         (DES_RADIUS(LL)+DES_RADIUS(I))
166                       IF(DIST_MAG > (VDW_INNER_CUTOFF+R_LM)**2) THEN
167                          DistApart = (SQRT(DIST_MAG)-R_LM)
168                          FORCE_COH = HAMAKER_CONSTANT * EQ_RADIUS /           &
169                             (12d0*DistApart**2) * (Asperities/(Asperities+    &
170                             EQ_RADIUS) + ONE/(ONE+Asperities/DistApart)**2 )
171                       ELSE
172                          FORCE_COH = 2d0 * PI * SURFACE_ENERGY * EQ_RADIUS *  &
173                            (Asperities/(Asperities+EQ_RADIUS) + ONE/          &
174                            (ONE+Asperities/VDW_INNER_CUTOFF)**2 )
175                       ENDIF
176                       FC_TMP(:) = DIST(:)*FORCE_COH/SQRT(DIST_MAG)
177                       TOW_TMP(:,:) = ZERO
178     
179     ! just for post-processing mag. of cohesive forces on each particle
180                       PostCohesive(LL) = PostCohesive(LL) + FORCE_COH / PMASS(LL)
181                    ENDIF
182                 ENDIF
183     
184                 IF (ENERGY_EQ) THEN
185                 ! Calculate conduction and radiation for thermodynamic neighbors
186                    IF(CALC_COND_DES(PIJK(LL,5))) THEN
187                       QQ_TMP = DES_CONDUCTION(LL, I, sqrt(DIST_MAG), PIJK(LL,5), PIJK(LL,4))
188     
189     !!$omp atomic
190                       Q_Source(LL) = Q_Source(LL) + QQ_TMP
191     
192     !!$omp atomic
193                       Q_Source(I) = Q_Source(I) - QQ_TMP
194                    ENDIF
195                 ENDIF
196     
197                 IF(DIST_MAG > (R_LM - SMALL_NUMBER)**2) THEN
198                    PFT_NEIGHBOR(:,CC) = 0.0
199                    CYCLE
200                 ENDIF
201     
202     #ifdef do_sap
203                    if (.not.is_pair(multisap%hashtable,ll,i)) then
204     
205                       print *,"SAP DIDNT FIND PAIR: ",ll,i
206                       print *,"PARTICLE (",ll,"):  ",des_pos_new(:,ll), " WITH RADIUS: ",des_radius(ll)
207                       print *,"PARTICLE (",i,"):  ",des_pos_new(:,i), " WITH RADIUS: ",des_radius(i)
208     
209                       print *,""
210                       print *," ******   ",sqrt(dot_product(des_pos_new(:,ll)-des_pos_new(:,i),des_pos_new(:,ll)-des_pos_new(:,i))),"     *********"
211                       print *,""
212     
213                       ! print *,"LLLLLLLLL ",boxhandle(ll)%list(:)
214                       ! print *,"IIIIIIIII ",boxhandle(i)%list(:)
215     
216                       do mm=1,size(boxhandle(ll)%list)
217                          if (boxhandle(ll)%list(mm)%sap_id < 0 ) cycle
218                          print *," PARTICLE ",ll," IS IN ",boxhandle(ll)%list(mm)%sap_id
219                          box_id = boxhandle(ll)%list(mm)%box_id
220     
221                          found = .false.
222                          do nn=1,size(boxhandle(i)%list)
223                             if (boxhandle(i)%list(nn)%sap_id .eq. boxhandle(ll)%list(mm)%sap_id) then
224                                ! print *," PARTICLE ",i," IS ALSO IN ",boxhandle(i)%list(nn)
225                                box_id2 = boxhandle(i)%list(nn)%box_id
226                                found = .true.
227                             endif
228                          enddo
229     
230                          if (.not.found) cycle
231     
232                          ! print *,"BOTH ",ll,i," ARE IN ",boxhandle(ll)%list(mm)
233     
234                       enddo
235     
236                       ERROR_STOP __LINE__
237                    endif
238     #endif
239     
240                 IF(DIST_MAG == 0) THEN
241                    WRITE(*,8550) LL, I
242                    ERROR_STOP "division by zero"
243      8550 FORMAT('distance between particles is zero:',2(2x,I10))
244                 ENDIF
245     
246                 DIST_MAG = SQRT(DIST_MAG)
247                 NORMAL(:)= DIST(:)/DIST_MAG
248     
249     ! Calcuate the normal overlap
250                 OVERLAP_N = R_LM-DIST_MAG
251                 IF(REPORT_EXCESS_OVERLAP) CALL PRINT_EXCESS_OVERLAP
252     
253     ! Calculate the components of translational relative velocity for a
254     ! contacting particle pair and the tangent to the plane of contact
255                 CALL CFRELVEL(LL, I, V_REL_TRANS_NORM, VREL_T,            &
256                    NORMAL(:), DIST_MAG)
257     
258                 phaseLL = PIJK(LL,5)
259                 phaseI = PIJK(I,5)
260     
261     ! Hertz spring-dashpot contact model
262                 IF (DES_COLL_MODEL_ENUM .EQ. HERTZIAN) THEN
263                    sqrt_overlap = SQRT(OVERLAP_N)
264                    KN_DES = hert_kn(phaseLL,phaseI)*sqrt_overlap
265                    KT_DES = hert_kt(phaseLL,phaseI)*sqrt_overlap
266                    sqrt_overlap = SQRT(sqrt_overlap)
267                    ETAN_DES = DES_ETAN(phaseLL,phaseI)*sqrt_overlap
268                    ETAT_DES = DES_ETAT(phaseLL,phaseI)*sqrt_overlap
269     
270     ! Linear spring-dashpot contact model
271                 ELSE
272                    KN_DES = KN
273                    KT_DES = KT
274                    ETAN_DES = DES_ETAN(phaseLL,phaseI)
275                    ETAT_DES = DES_ETAT(phaseLL,phaseI)
276                 ENDIF
277     
278     ! Calculate the normal contact force
279                 FN(:) =  -(KN_DES * OVERLAP_N * NORMAL(:) + &
280                    ETAN_DES * V_REL_TRANS_NORM * NORMAL(:))
281     
282     ! Calcuate the tangential overlap
283                 OVERLAP_T(:) = DTSOLID*VREL_T(:) + PFT_NEIGHBOR(:,CC)
284                 MAG_OVERLAP_T = sqrt(dot_product(OVERLAP_T,OVERLAP_T))
285     
286     ! Calculate the tangential contact force.
287                 IF(MAG_OVERLAP_T > 0.0) THEN
288     ! Tangential froce from spring.
289                    FTMD = KT_DES*MAG_OVERLAP_T
290     ! Max force before the on set of frictional slip.
291                    FNMD = MEW*sqrt(dot_product(FN,FN))
292     ! Direction of tangential force.
293                    TANGENT = OVERLAP_T/MAG_OVERLAP_T
294     ! Frictional slip
295                    IF(FTMD > FNMD) THEN
296                       FT = -FNMD * TANGENT
297                       OVERLAP_T = (FNMD/KT_DES) * TANGENT
298                    ELSE
299                       FT = -FTMD * TANGENT
300                    ENDIF
301                 ELSE
302                    FT = 0.0
303                 ENDIF
304     
305     ! Add in the tangential dashpot damping force
306                 FT = FT - ETAT_DES * VREL_T(:)
307     
308     ! Save tangential displacement history
309                 PFT_NEIGHBOR(:,CC) = OVERLAP_T(:)
310     
311     ! calculate the distance from the particles' centers to the contact point,
312     ! which is taken as the radical line
313     ! dist_ci+dist_cl=dist_li; dist_ci^2+a^2=ri^2;  dist_cl^2+a^2=rl^2
314                 DIST_CL = DIST_MAG/2.d0 + (DES_RADIUS(LL)**2 - &
315                    DES_RADIUS(I)**2)/(2.d0*DIST_MAG)
316     
317                 DIST_CI = DIST_MAG - DIST_CL
318     
319                 TOW_force(:) = DES_CROSSPRDCT(NORMAL(:), FT(:))
320                 TOW_TMP(:,1) = DIST_CL*TOW_force(:)
321                 TOW_TMP(:,2) = DIST_CI*TOW_force(:)
322     
323     ! Calculate the total force FC of a collision pair
324     ! total contact force ( FC_TMP may already include cohesive force)
325                 FC_TMP(:) = FC_TMP(:) + FN(:) + FT(:)
326     
327                 FC(LL,:) = FC(LL,:) + FC_TMP(:)
328     
329     !!$omp atomic
330                 FC(I,1) = FC(I,1) - FC_TMP(1)
331     !!$omp atomic
332                 FC(I,2) = FC(I,2) - FC_TMP(2)
333     !!$omp atomic
334                 FC(I,3) = FC(I,3) - FC_TMP(3)
335     
336     ! for each particle the signs of norm and ft both flip, so add the same torque
337                 TOW(LL,:) = TOW(LL,:) + TOW_TMP(:,1)
338     
339     !!$omp atomic
340                 TOW(I,1)  = TOW(I,1)  + TOW_TMP(1,2)
341     !!$omp atomic
342                 TOW(I,2)  = TOW(I,2)  + TOW_TMP(2,2)
343     !!$omp atomic
344                 TOW(I,3)  = TOW(I,3)  + TOW_TMP(3,2)
345     
346              ENDDO
347           ENDDO
348     !!$omp end do
349     
350     !!$omp end parallel
351     
352     ! just for post-processing mag. of cohesive forces on each particle
353           IF(USE_COHESION .AND. VAN_DER_WAALS .AND. GRAV_MAG > ZERO) THEN
354              PostCohesive(:) = PostCohesive(:)/GRAV_MAG
355           ENDIF
356     
357           RETURN
358     
359           contains
360     
361             include 'functions.inc'
362     
363     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
364     !                                                                      !
365     !  Subroutine: print_excess_overlap                                    !
366     !                                                                      !
367     !  Purpose: Print overlap warning messages.                            !
368     !                                                                      !
369     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
370           SUBROUTINE PRINT_EXCESS_OVERLAP
371     
372           use error_manager
373     
374           IF(OVERLAP_N > flag_overlap*DES_RADIUS(LL) .OR.                  &
375              OVERLAP_N > flag_overlap*DES_RADIUS(I)) THEN
376     
377              WRITE(ERR_MSG,1000) trim(iVAL(LL)), trim(iVAL(I)), S_TIME,    &
378                 DES_RADIUS(LL), DES_RADIUS(I), OVERLAP_N
379     
380              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
381           ENDIF
382     
383      1000 FORMAT('WARNING: Excessive overplay detected between ',          &
384              'particles ',A,' and ',/A,' at time ',g11.4,'.',/             &
385              'RADII:  ',g11.4,' and ',g11.4,4x,'OVERLAP: ',g11.4)
386     
387           END SUBROUTINE PRINT_EXCESS_OVERLAP
388     
389         END SUBROUTINE CALC_FORCE_DEM
390