File: /nfs/home/0/users/jenkins/mfix.git/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           SUBROUTINE CALC_FORCE_DEM
13     
14     !---------------------------------------------------------------------//
15           USE calc_collision_wall
16           USE constant, ONLY: Pi
17           USE des_thermo
18           USE des_thermo_cond
19           USE discretelement
20           USE physprop, ONLY: K_s0
21           USE run
22     
23           IMPLICIT NONE
24     
25     ! Local variables
26     !---------------------------------------------------------------------//
27     ! percent of particle radius when excess overlap will be flagged
28           DOUBLE PRECISION, PARAMETER :: flag_overlap = 0.20d0
29     ! particle no. indices
30           INTEGER :: I, LL, CC
31     ! the overlap occuring between particle-particle or particle-wall
32     ! collision in the normal direction
33           DOUBLE PRECISION :: OVERLAP_N
34     ! square root of the overlap
35           DOUBLE PRECISION :: SQRT_OVERLAP
36     ! distance vector between two particle centers or between a particle
37     ! center and wall when the two surfaces are just at contact (i.e. no
38     ! overlap)
39           DOUBLE PRECISION :: R_LM,DIST_CI,DIST_CL
40     ! the normal and tangential components of the translational relative
41     ! velocity
42           DOUBLE PRECISION :: V_REL_TRANS_NORM
43     ! distance vector between two particle centers or between a particle
44     ! center and wall at current and previous time steps
45           DOUBLE PRECISION :: DIST(3), DIST_NORM(3), DIST_MAG
46     ! tangent to the plane of contact at current time step
47           DOUBLE PRECISION :: V_REL_TANG(3)
48     ! normal and tangential forces
49           DOUBLE PRECISION :: FN(3), FT(3)
50           DOUBLE PRECISION :: FNS1(3), FNS2(3)
51           DOUBLE PRECISION :: FTS1(3), FTS2(3)
52     ! temporary storage of tangential DISPLACEMENT
53           DOUBLE PRECISION :: PFT_TMP(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     ! set to T when a sliding contact occurs
71           LOGICAL :: PARTICLE_SLIDE
72     
73           LOGICAL, PARAMETER :: report_excess_overlap = .FALSE.
74     
75     !-----------------------------------------------
76     
77     ! Initialize cohesive forces
78           IF(USE_COHESION) PostCohesive(:) = ZERO
79     
80           CALL CALC_DEM_FORCE_WITH_WALL_STL
81     
82     
83     ! Check particle LL neighbour contacts
84     !---------------------------------------------------------------------//
85     
86     !$omp parallel default(none) private(cc,ll,i,dist,r_lm,             &
87     !$omp    overlap_n,v_rel_tang,v_rel_trans_norm,sqrt_overlap,           &
88     !$omp    kn_des,kt_des,hert_kn,hert_kt,phasell,phasei,etan_des,        &
89     !$omp    etat_des,fns1,fns2,fts1,fts2,pft_tmp,fn,ft,particle_slide,    &
90     !$omp    eq_radius,distapart,force_coh,k_s0,dist_mag,dist_norm,        &
91     !$omp    dist_cl, dist_ci, fc_tmp, tow_tmp, tow_force, qq_tmp)         &
92     !$omp    shared(pairs,pair_num,des_pos_new,des_radius,                 &
93     !$omp    des_coll_model_enum,kn,kt,pv_pair,pft_pair,pfn_pair,pijk,     &
94     !$omp    des_etan,des_etat,mew,use_cohesion, calc_cond_des,            &
95     !$omp    van_der_waals,vdw_outer_cutoff,vdw_inner_cutoff,              &
96     !$omp    hamaker_constant,asperities,surface_energy, pea,              &
97     !$omp    tow, fc, energy_eq, grav_mag, postcohesive, pmass, q_source)
98     
99     !$omp do
100           DO CC = 1, PAIR_NUM
101              LL = PAIRS(1,CC)
102              I  = PAIRS(2,CC)
103     
104              IF(.NOT.PEA(LL,1)) CYCLE
105              IF(.NOT.PEA(I, 1)) CYCLE
106     
107              R_LM = DES_RADIUS(LL) + DES_RADIUS(I)
108              DIST(:) = DES_POS_NEW(:,I) - DES_POS_NEW(:,LL)
109              DIST_MAG = dot_product(DIST,DIST)
110     
111              FC_TMP(:) = ZERO
112     
113     ! Compute particle-particle VDW cohesive short-range forces
114              IF(USE_COHESION .AND. VAN_DER_WAALS) THEN
115                 IF(DIST_MAG < (R_LM+VDW_OUTER_CUTOFF)**2) THEN
116                    EQ_RADIUS = 2d0 * DES_RADIUS(LL)*DES_RADIUS(I) /        &
117                         (DES_RADIUS(LL)+DES_RADIUS(I))
118                    IF(DIST_MAG > (VDW_INNER_CUTOFF+R_LM)**2) THEN
119                       DistApart = (SQRT(DIST_MAG)-R_LM)
120                       FORCE_COH = HAMAKER_CONSTANT * EQ_RADIUS /           &
121                          (12d0*DistApart**2) * (Asperities/(Asperities+    &
122                          EQ_RADIUS) + ONE/(ONE+Asperities/DistApart)**2 )
123                    ELSE
124                       FORCE_COH = 2d0 * PI * SURFACE_ENERGY * EQ_RADIUS *  &
125                         (Asperities/(Asperities+EQ_RADIUS) + ONE/          &
126                         (ONE+Asperities/VDW_INNER_CUTOFF)**2 )
127                    ENDIF
128                    FC_TMP(:) = DIST(:)*FORCE_COH/SQRT(DIST_MAG)
129                    TOW_TMP(:,:) = ZERO
130     
131     ! just for post-processing mag. of cohesive forces on each particle
132     
133                    PostCohesive(LL) = PostCohesive(LL)
134                    if(GRAV_MAG > ZERO .AND. PEA(LL,1)) THEN
135                       FORCE_COH = SQRT(dot_product(FC_TMP(:),FC_TMP(:))) / (PMASS(LL)*GRAV_MAG)
136     
137                       !$omp atomic
138                       PostCohesive(LL) = PostCohesive(LL) + FORCE_COH
139                    ENDIF
140                 ENDIF
141              ENDIF
142     
143              IF (ENERGY_EQ) THEN
144                 ! Calculate conduction and radiation for thermodynamic neighbors
145                 IF(CALC_COND_DES(PIJK(LL,5))) THEN
146                    QQ_TMP = DES_CONDUCTION(LL, I, sqrt(DIST_MAG), PIJK(LL,5), PIJK(LL,4))
147     
148                    !$omp atomic
149                    Q_Source(LL) = Q_Source(LL) + QQ_TMP
150     
151                    !$omp atomic
152                    Q_Source(I) = Q_Source(I) - QQ_TMP
153                 ENDIF
154              ENDIF
155     
156              IF(DIST_MAG > (R_LM + SMALL_NUMBER)**2) THEN
157                 PV_PAIR(CC) = .false.
158                 PFT_PAIR(:,CC) = 0.0
159                 PFN_PAIR(:,CC) = 0.0
160                 CYCLE
161              ENDIF
162     
163              IF(DIST_MAG == 0) THEN
164                 WRITE(*,8550) LL, I
165                 STOP "division by zero"
166      8550 FORMAT('distance between particles is zero:',2(2x,I10))
167              ENDIF
168              DIST_MAG = SQRT(DIST_MAG)
169              DIST_NORM(:)= DIST(:)/DIST_MAG
170     
171     ! Overlap calculation changed from history based to current position
172              OVERLAP_N = R_LM-DIST_MAG
173     
174              IF (report_excess_overlap) call print_excess_overlap
175     
176     ! Calculate the components of translational relative velocity for a
177     ! contacting particle pair and the tangent to the plane of contact
178              CALL CFRELVEL(LL, I, V_REL_TRANS_NORM, V_REL_TANG,            &
179                 DIST_NORM(:), DIST_MAG)
180     
181              phaseLL = PIJK(LL,5)
182              phaseI = PIJK(I,5)
183     
184     ! Hertz spring-dashpot contact model
185              IF (DES_COLL_MODEL_ENUM .EQ. HERTZIAN) THEN
186                 sqrt_overlap = SQRT(OVERLAP_N)
187                 KN_DES = hert_kn(phaseLL,phaseI)*sqrt_overlap
188                 KT_DES = hert_kt(phaseLL,phaseI)*sqrt_overlap
189                 sqrt_overlap = SQRT(sqrt_overlap)
190                 ETAN_DES = DES_ETAN(phaseLL,phaseI)*sqrt_overlap
191                 ETAT_DES = DES_ETAT(phaseLL,phaseI)*sqrt_overlap
192     
193     ! Linear spring-dashpot contact model
194              ELSE
195                 KN_DES = KN
196                 KT_DES = KT
197                 ETAN_DES = DES_ETAN(phaseLL,phaseI)
198                 ETAT_DES = DES_ETAT(phaseLL,phaseI)
199              ENDIF
200     
201     ! Calculate the normal contact force
202              FNS1(:) = -KN_DES * OVERLAP_N * DIST_NORM(:)
203              FNS2(:) = -ETAN_DES * V_REL_TRANS_NORM*DIST_NORM(:)
204              FN(:) = FNS1(:) + FNS2(:)
205     
206              call calc_tangential_displacement(pft_tmp(:),DIST_NORM(:), &
207                 pfn_pair(:,cc),pft_pair(:,cc),overlap_n,                   &
208                 v_rel_trans_norm,v_rel_tang(:),PV_PAIR(CC))
209              PV_PAIR(CC) = .true.
210     
211     ! Calculate the tangential contact force
212              FTS1(:) = -KT_DES * PFT_TMP(:)
213              FTS2(:) = -ETAT_DES * V_REL_TANG
214              FT(:) = FTS1(:) + FTS2(:)
215     
216     ! Check for Coulombs friction law and limit the maximum value of the
217     ! tangential force on a particle in contact with another particle/wall
218              PARTICLE_SLIDE = .FALSE.
219              CALL CFSLIDE(V_REL_TANG(:), PARTICLE_SLIDE, MEW,              &
220                  FT(:), FN(:))
221     ! calculate the distance from the particles' centers to the contact point,
222     ! which is taken as the radical line
223     ! dist_ci+dist_cl=dist_li; dist_ci^2+a^2=ri^2;  dist_cl^2+a^2=rl^2
224              DIST_CL = DIST_MAG/2.d0 + (DES_RADIUS(LL)**2 - &
225                   DES_RADIUS(I)**2)/(2.d0*DIST_MAG)
226     
227              DIST_CI = DIST_MAG - DIST_CL
228     
229              TOW_force(:) = DES_CROSSPRDCT(DIST_NORM(:), FT(:))
230              TOW_TMP(:,1) = DIST_CL*TOW_force(:)
231              TOW_TMP(:,2) = DIST_CI*TOW_force(:)
232     
233     ! Calculate the total force FC of a collision pair
234     ! total contact force ( FC_TMP may already include cohesive force)
235              FC_TMP(:) = FC_TMP(:) + FN(:) + FT(:)
236                 !$omp atomic
237                 FC(1,LL) = FC(1,LL) + FC_TMP(1)
238                 !$omp atomic
239                 FC(2,LL) = FC(2,LL) + FC_TMP(2)
240                 !$omp atomic
241                 FC(3,LL) = FC(3,LL) + FC_TMP(3)
242     
243                 I  = PAIRS(2,CC)
244                 !$omp atomic
245                 FC(1,I) = FC(1,I) - FC_TMP(1)
246                 !$omp atomic
247                 FC(2,I) = FC(2,I) - FC_TMP(2)
248                 !$omp atomic
249                 FC(3,I) = FC(3,I) - FC_TMP(3)
250     
251     
252     ! for each particle the signs of norm and ft both flip, so add the same torque
253                 !$omp atomic
254                 TOW(1,LL) = TOW(1,LL) + TOW_TMP(1,1)
255                 !$omp atomic
256                 TOW(2,LL) = TOW(2,LL) + TOW_TMP(2,1)
257                 !$omp atomic
258                 TOW(3,LL) = TOW(3,LL) + TOW_TMP(3,1)
259     
260                 !$omp atomic
261                 TOW(1,I)  = TOW(1,I)  + TOW_TMP(1,2)
262                 !$omp atomic
263                 TOW(2,I)  = TOW(2,I)  + TOW_TMP(2,2)
264                 !$omp atomic
265                 TOW(3,I)  = TOW(3,I)  + TOW_TMP(3,2)
266     
267     
268     ! Save tangential displacement history with Coulomb's law correction
269              IF (PARTICLE_SLIDE) THEN
270     ! Since FT might be corrected during the call to cfslide, the tangential
271     ! displacement history needs to be changed accordingly
272                 PFT_PAIR(:,CC) = -( FT(:) - FTS2(:) ) / KT_DES
273              ELSE
274                 PFT_PAIR(:,CC) = PFT_TMP(:)
275              ENDIF
276     
277           ENDDO
278     !$omp end do
279     
280     !$omp end parallel
281     
282           RETURN
283     
284           contains
285     
286     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
287     !                                                                      !
288     !  Subroutine: print_excess_overlap                                    !
289     !                                                                      !
290     !  Purpose: Print overlap warning messages.                            !
291     !                                                                      !
292     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
293           SUBROUTINE PRINT_EXCESS_OVERLAP
294     
295           use error_manager
296     
297           IF(OVERLAP_N > flag_overlap*DES_RADIUS(LL) .OR.                  &
298              OVERLAP_N > flag_overlap*DES_RADIUS(I)) THEN
299     
300              WRITE(ERR_MSG,1000) trim(iVAL(LL)), trim(iVAL(I)), S_TIME,    &
301                 DES_RADIUS(LL), DES_RADIUS(I), OVERLAP_N
302     
303              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
304           ENDIF
305     
306      1000 FORMAT('WARNING: Excessive overplay detected between ',          &
307              'particles ',A,' and ',/A,' at time ',g11.4,'.',/             &
308              'RADII:  ',g11.4,' and ',g11.4,4x,'OVERLAP: ',g11.4)
309     
310           END SUBROUTINE PRINT_EXCESS_OVERLAP
311     
312     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
313     !                                                                      !
314     !  Subroutine: CALC_TANGENTIAL_DISPLACEMENT                            !
315     !                                                                      !
316     !  Purpose: Calculate the tangential displacement which is the         !
317     !  integration of tangential relative velocity with respect to         !
318     !  contact time.                                                       !
319     !                                                                      !
320     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
321           SUBROUTINE CALC_TANGENTIAL_DISPLACEMENT(PFT, NORM, NORM_OLD,     &
322              SIGMAT_OLD, OVERLAP_NORM, RELVEL_TANG_NORM, RELVEL_TANG,      &
323              ALREADY_COLLIDING)
324     
325           IMPLICIT NONE
326     
327     ! tangential displacement
328           DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: PFT
329     
330     ! local variables for the accumulated tangential displacement that occurs
331     ! for the particle-particle or particle-wall collision (current time
332     ! step and previous time step)
333           DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: SIGMAT_OLD
334           LOGICAL, INTENT(IN) :: ALREADY_COLLIDING
335     
336     ! the overlap occuring between particle-particle or particle-wall
337     ! collision in the tangential direction
338           DOUBLE PRECISION :: OVERLAP_T(3)
339     
340           DOUBLE PRECISION, DIMENSION(3) :: SIGMAT
341     
342     ! time elapsed to travel the calculated normal overlap given the
343     ! normal relative velocity
344           DOUBLE PRECISION :: DTSOLID_TMP
345     
346     ! unit normal vector along the line of contact between contacting
347     ! particles or particle-wall at previous time step
348           DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: NORM, NORM_OLD
349           DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: RELVEL_TANG
350           DOUBLE PRECISION, INTENT(IN) :: OVERLAP_NORM, RELVEL_TANG_NORM
351     
352     ! variables for tangential displacement calculation:
353     ! unit vector for axis of rotation and its magnitude
354           DOUBLE PRECISION :: TMP_AX(3), TMP_MAG
355     
356     ! tangent to the plane of contact at current and previous time step
357     ! (used for 2D calculations)
358           DOUBLE PRECISION :: TANG_OLD(3),TANG_NEW(3)
359     
360           IF(already_colliding) THEN
361              OVERLAP_T(:) = DTSOLID*relvel_tang(:)
362           ELSE
363              DTSOLID_TMP = OVERLAP_NORM/MAX(relvel_tang_norm,SMALL_NUMBER)
364              OVERLAP_T(:) = MIN(DTSOLID,DTSOLID_TMP)*relvel_tang(:)
365           ENDIF
366     
367           IF(USE_VDH_DEM_MODEL) then
368     ! Calculate the tangential displacement which is integration of
369     ! tangential relative velocity with respect to contact time.
370     ! Correction in the tangential direction is imposed
371     
372     ! New procedure: van der Hoef et al. (2006)
373     
374     ! calculate the unit vector for axis of rotation
375              tmp_ax = des_crossprdct(norm_old,norm)
376              tmp_mag=dot_product(tmp_ax,tmp_ax)
377              if(tmp_mag .gt. zero)then
378                 tmp_ax(:)=tmp_ax(:)/sqrt(tmp_mag)
379     ! get the old tangential direction unit vector
380                 tang_old = des_crossprdct(tmp_ax,norm_old)
381     ! get the new tangential direction unit vector due to rotation
382                 tang_new = des_crossprdct(tmp_ax,norm)
383                 sigmat(:)=dot_product(sigmat_old,tmp_ax)*tmp_ax(:) &
384                      + dot_product(sigmat_old,tang_old)*tang_new(:)
385                 sigmat(:)=sigmat(:)+overlap_t(:)
386              else
387                 sigmat(:)=sigmat_old(:)+overlap_t(:)
388              endif
389     
390     ! Save the old normal direction
391              PFT(:)   = SIGMAT(:)
392           ELSE
393              ! Old procedure
394              sigmat = sigmat_old + OVERLAP_T(:)
395              pft(:) = sigmat - dot_product(sigmat,NORM)*NORM(:)
396           ENDIF
397     
398         END SUBROUTINE CALC_TANGENTIAL_DISPLACEMENT
399     
400         END SUBROUTINE CALC_FORCE_DEM
401