File: RELATIVE:/../../../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 run
21     
22           IMPLICIT NONE
23     
24     ! Local variables
25     !---------------------------------------------------------------------//
26     ! percent of particle radius when excess overlap will be flagged
27           DOUBLE PRECISION, PARAMETER :: flag_overlap = 0.20d0
28     ! particle no. indices
29           INTEGER :: I, LL, cc, CC_START, CC_END
30     ! the overlap occuring between particle-particle or particle-wall
31     ! collision in the normal direction
32           DOUBLE PRECISION :: OVERLAP_N, OVERLAP_T(3)
33     ! square root of the overlap
34           DOUBLE PRECISION :: SQRT_OVERLAP
35     ! distance vector between two particle centers or between a particle
36     ! center and wall when the two surfaces are just at contact (i.e. no
37     ! overlap)
38           DOUBLE PRECISION :: R_LM,DIST_CI,DIST_CL
39     ! the normal and tangential components of the translational relative
40     ! velocity
41           DOUBLE PRECISION :: V_REL_TRANS_NORM, rad
42     ! distance vector between two particle centers or between a particle
43     ! center and wall at current and previous time steps
44           DOUBLE PRECISION :: DIST(3), NORMAL(3), DIST_MAG, POS(3)
45     ! tangent to the plane of contact at current time step
46           DOUBLE PRECISION :: VREL_T(3)
47     ! normal and tangential forces
48           DOUBLE PRECISION :: FN(3), FT(3)
49     ! temporary storage of force
50           DOUBLE PRECISION :: FC_TMP(3)
51     ! temporary storage of force for torque
52           DOUBLE PRECISION :: TOW_FORCE(3)
53     ! temporary storage of torque
54           DOUBLE PRECISION :: TOW_TMP(3,2)
55     ! temporary storage of conduction/radiation
56           DOUBLE PRECISION :: QQ_TMP
57     
58     ! store solids phase index of particle (i.e. pijk(np,5))
59           INTEGER :: PHASEI, PHASELL
60     ! local values used spring constants and damping coefficients
61           DOUBLE PRECISION :: ETAN_DES, ETAT_DES
62           DOUBLE PRECISION :: KN_DES, KT_DES
63     ! local values used for calculating cohesive forces
64           DOUBLE PRECISION :: FORCE_COH, EQ_RADIUS, DistApart
65     
66           LOGICAL, PARAMETER :: report_excess_overlap = .FALSE.
67     
68           DOUBLE PRECISION :: FNMD, FTMD, MAG_OVERLAP_T, TANGENT(3)
69     
70     !-----------------------------------------------
71     
72     ! Initialize cohesive forces
73           IF(USE_COHESION) PostCohesive(:) = ZERO
74     
75           CALL CALC_DEM_FORCE_WITH_WALL_STL
76     
77     ! Check particle LL neighbor contacts
78     !---------------------------------------------------------------------//
79     
80     !$omp parallel default(none) private(pos,rad,cc,cc_start,cc_end,ll,i,  &
81     !$omp    overlap_n,vrel_t,v_rel_trans_norm,sqrt_overlap,dist,r_lm,     &
82     !$omp    kn_des,kt_des,hert_kn,hert_kt,phasell,phasei,etan_des,        &
83     !$omp    etat_des,fn,ft,overlap_t,tangent,mag_overlap_t,               &
84     !$omp    eq_radius,distapart,force_coh,dist_mag,NORMAL,ftmd,fnmd,      &
85     !$omp    dist_cl, dist_ci, fc_tmp, tow_tmp, tow_force, qq_tmp)         &
86     !$omp shared(max_pip,neighbors,neighbor_index,des_pos_new,des_radius,  &
87     !$omp    des_coll_model_enum,kn,kt,pft_neighbor,pijk,neigh_max,        &
88     !$omp    des_etan,des_etat,mew,use_cohesion, calc_cond_des, dtsolid,   &
89     !$omp    van_der_waals,vdw_outer_cutoff,vdw_inner_cutoff,              &
90     !$omp    hamaker_constant,asperities,surface_energy,                   &
91     !$omp    tow, fc, energy_eq, grav_mag, postcohesive, pmass, q_source)
92     
93     !$omp do
94     
95           DO LL = 1, MAX_PIP
96              IF(IS_NONEXISTENT(LL)) CYCLE
97              pos = DES_POS_NEW(:,LL)
98              rad = DES_RADIUS(LL)
99     
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                 IF(IS_NONEXISTENT(I)) CYCLE
107     
108                 R_LM = rad + DES_RADIUS(I)
109                 DIST(:) = DES_POS_NEW(:,I) - POS(:)
110                 DIST_MAG = dot_product(DIST,DIST)
111     
112                 FC_TMP(:) = ZERO
113     
114     ! Compute particle-particle VDW cohesive short-range forces
115                 IF(USE_COHESION .AND. VAN_DER_WAALS) THEN
116                    IF(DIST_MAG < (R_LM+VDW_OUTER_CUTOFF)**2) THEN
117                       EQ_RADIUS = 2d0 * DES_RADIUS(LL)*DES_RADIUS(I) /     &
118                         (DES_RADIUS(LL)+DES_RADIUS(I))
119                       IF(DIST_MAG > (VDW_INNER_CUTOFF+R_LM)**2) THEN
120                          DistApart = (SQRT(DIST_MAG)-R_LM)
121                          FORCE_COH = HAMAKER_CONSTANT * EQ_RADIUS /           &
122                             (12d0*DistApart**2) * (Asperities/(Asperities+    &
123                             EQ_RADIUS) + ONE/(ONE+Asperities/DistApart)**2 )
124                       ELSE
125                          FORCE_COH = 2d0 * PI * SURFACE_ENERGY * EQ_RADIUS *  &
126                            (Asperities/(Asperities+EQ_RADIUS) + ONE/          &
127                            (ONE+Asperities/VDW_INNER_CUTOFF)**2 )
128                       ENDIF
129                       FC_TMP(:) = DIST(:)*FORCE_COH/SQRT(DIST_MAG)
130                       TOW_TMP(:,:) = ZERO
131     
132     ! just for post-processing mag. of cohesive forces on each particle
133                       PostCohesive(LL) = PostCohesive(LL) + FORCE_COH / PMASS(LL)
134                    ENDIF
135                 ENDIF
136     
137                 IF (ENERGY_EQ) THEN
138                 ! Calculate conduction and radiation for thermodynamic neighbors
139                    IF(CALC_COND_DES(PIJK(LL,5))) THEN
140                       QQ_TMP = DES_CONDUCTION(LL, I, sqrt(DIST_MAG), PIJK(LL,5), PIJK(LL,4))
141     
142                       !$omp atomic
143                       Q_Source(LL) = Q_Source(LL) + QQ_TMP
144     
145                       !$omp atomic
146                       Q_Source(I) = Q_Source(I) - QQ_TMP
147                    ENDIF
148                 ENDIF
149     
150                 IF(DIST_MAG > (R_LM + SMALL_NUMBER)**2) THEN
151                    PFT_NEIGHBOR(:,CC) = 0.0
152                    CYCLE
153                 ENDIF
154     
155                 IF(DIST_MAG == 0) THEN
156                    WRITE(*,8550) LL, I
157                    STOP "division by zero"
158      8550 FORMAT('distance between particles is zero:',2(2x,I10))
159                 ENDIF
160     
161                 DIST_MAG = SQRT(DIST_MAG)
162                 NORMAL(:)= DIST(:)/DIST_MAG
163     
164     ! Calcuate the normal overlap
165                 OVERLAP_N = R_LM-DIST_MAG
166                 IF(REPORT_EXCESS_OVERLAP) CALL PRINT_EXCESS_OVERLAP
167     
168     ! Calculate the components of translational relative velocity for a
169     ! contacting particle pair and the tangent to the plane of contact
170                 CALL CFRELVEL(LL, I, V_REL_TRANS_NORM, VREL_T,            &
171                    NORMAL(:), DIST_MAG)
172     
173                 phaseLL = PIJK(LL,5)
174                 phaseI = PIJK(I,5)
175     
176     ! Hertz spring-dashpot contact model
177                 IF (DES_COLL_MODEL_ENUM .EQ. HERTZIAN) THEN
178                    sqrt_overlap = SQRT(OVERLAP_N)
179                    KN_DES = hert_kn(phaseLL,phaseI)*sqrt_overlap
180                    KT_DES = hert_kt(phaseLL,phaseI)*sqrt_overlap
181                    sqrt_overlap = SQRT(sqrt_overlap)
182                    ETAN_DES = DES_ETAN(phaseLL,phaseI)*sqrt_overlap
183                    ETAT_DES = DES_ETAT(phaseLL,phaseI)*sqrt_overlap
184     
185     ! Linear spring-dashpot contact model
186                 ELSE
187                    KN_DES = KN
188                    KT_DES = KT
189                    ETAN_DES = DES_ETAN(phaseLL,phaseI)
190                    ETAT_DES = DES_ETAT(phaseLL,phaseI)
191                 ENDIF
192     
193     ! Calculate the normal contact force
194                 FN(:) =  -(KN_DES * OVERLAP_N * NORMAL(:) + &
195                    ETAN_DES * V_REL_TRANS_NORM * NORMAL(:))
196     
197     ! Calcuate the tangential overlap
198                 OVERLAP_T(:) = DTSOLID*VREL_T(:) + PFT_NEIGHBOR(:,CC)
199                 MAG_OVERLAP_T = sqrt(dot_product(OVERLAP_T,OVERLAP_T))
200     
201     ! Calculate the tangential contact force.
202                 IF(MAG_OVERLAP_T > 0.0) THEN
203     ! Tangential froce from spring.
204                    FTMD = KT_DES*MAG_OVERLAP_T
205     ! Max force before the on set of frictional slip.
206                    FNMD = MEW*sqrt(dot_product(FN,FN))
207     ! Direction of tangential force.
208                    TANGENT = OVERLAP_T/MAG_OVERLAP_T
209     ! Frictional slip
210                    IF(FTMD > FNMD) THEN
211                       FT = -FNMD * TANGENT
212                       OVERLAP_T = (FNMD/KT_DES) * TANGENT
213                    ELSE
214                       FT = -FTMD * TANGENT
215                    ENDIF
216                 ELSE
217                    FT = 0.0
218                 ENDIF
219     
220     ! Add in the tangential dashpot damping force
221                 FT = FT - ETAT_DES * VREL_T(:)
222     
223     ! Save tangential displacement history
224                 PFT_NEIGHBOR(:,CC) = OVERLAP_T(:)
225     
226     ! calculate the distance from the particles' centers to the contact point,
227     ! which is taken as the radical line
228     ! dist_ci+dist_cl=dist_li; dist_ci^2+a^2=ri^2;  dist_cl^2+a^2=rl^2
229                 DIST_CL = DIST_MAG/2.d0 + (DES_RADIUS(LL)**2 - &
230                    DES_RADIUS(I)**2)/(2.d0*DIST_MAG)
231     
232                 DIST_CI = DIST_MAG - DIST_CL
233     
234                 TOW_force(:) = DES_CROSSPRDCT(NORMAL(:), FT(:))
235                 TOW_TMP(:,1) = DIST_CL*TOW_force(:)
236                 TOW_TMP(:,2) = DIST_CI*TOW_force(:)
237     
238     ! Calculate the total force FC of a collision pair
239     ! total contact force ( FC_TMP may already include cohesive force)
240                 FC_TMP(:) = FC_TMP(:) + FN(:) + FT(:)
241     
242                 FC(:,LL) = FC(:,LL) + FC_TMP(:)
243     
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     ! for each particle the signs of norm and ft both flip, so add the same torque
252                 TOW(:,LL) = TOW(:,LL) + TOW_TMP(:,1)
253     
254                 !$omp atomic
255                 TOW(1,I)  = TOW(1,I)  + TOW_TMP(1,2)
256                 !$omp atomic
257                 TOW(2,I)  = TOW(2,I)  + TOW_TMP(2,2)
258                 !$omp atomic
259                 TOW(3,I)  = TOW(3,I)  + TOW_TMP(3,2)
260     
261              ENDDO
262           ENDDO
263     !$omp end do
264     
265     !$omp end parallel
266     
267     ! just for post-processing mag. of cohesive forces on each particle
268           IF(USE_COHESION .AND. VAN_DER_WAALS .AND. GRAV_MAG > ZERO) THEN
269              PostCohesive(:) = PostCohesive(:)/GRAV_MAG
270           ENDIF
271     
272           RETURN
273     
274           contains
275     
276             include 'functions.inc'
277     
278     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
279     !                                                                      !
280     !  Subroutine: print_excess_overlap                                    !
281     !                                                                      !
282     !  Purpose: Print overlap warning messages.                            !
283     !                                                                      !
284     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
285           SUBROUTINE PRINT_EXCESS_OVERLAP
286     
287           use error_manager
288     
289           IF(OVERLAP_N > flag_overlap*DES_RADIUS(LL) .OR.                  &
290              OVERLAP_N > flag_overlap*DES_RADIUS(I)) THEN
291     
292              WRITE(ERR_MSG,1000) trim(iVAL(LL)), trim(iVAL(I)), S_TIME,    &
293                 DES_RADIUS(LL), DES_RADIUS(I), OVERLAP_N
294     
295              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
296           ENDIF
297     
298      1000 FORMAT('WARNING: Excessive overplay detected between ',          &
299              'particles ',A,' and ',/A,' at time ',g11.4,'.',/             &
300              'RADII:  ',g11.4,' and ',g11.4,4x,'OVERLAP: ',g11.4)
301     
302           END SUBROUTINE PRINT_EXCESS_OVERLAP
303     
304         END SUBROUTINE CALC_FORCE_DEM
305