File: /nfs/home/0/users/jenkins/mfix.git/model/des/calc_force_dem.f
1
2
3
4
5
6
7
8
9
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
26
27
28 DOUBLE PRECISION, PARAMETER :: flag_overlap = 0.20d0
29
30 INTEGER :: I, LL, CC
31
32
33 DOUBLE PRECISION :: OVERLAP_N
34
35 DOUBLE PRECISION :: SQRT_OVERLAP
36
37
38
39 DOUBLE PRECISION :: R_LM,DIST_CI,DIST_CL
40
41
42 DOUBLE PRECISION :: V_REL_TRANS_NORM
43
44
45 DOUBLE PRECISION :: DIST(3), DIST_NORM(3), DIST_MAG
46
47 DOUBLE PRECISION :: V_REL_TANG(3)
48
49 DOUBLE PRECISION :: FN(3), FT(3)
50 DOUBLE PRECISION :: FNS1(3), FNS2(3)
51 DOUBLE PRECISION :: FTS1(3), FTS2(3)
52
53 DOUBLE PRECISION :: PFT_TMP(3)
54
55 DOUBLE PRECISION :: FC_TMP(3)
56
57 DOUBLE PRECISION :: TOW_FORCE(3)
58
59 DOUBLE PRECISION :: TOW_TMP(3,2)
60
61 DOUBLE PRECISION :: QQ_TMP
62
63
64 INTEGER :: PHASEI, PHASELL
65
66 DOUBLE PRECISION :: ETAN_DES, ETAT_DES
67 DOUBLE PRECISION :: KN_DES, KT_DES
68
69 DOUBLE PRECISION :: FORCE_COH, EQ_RADIUS, DistApart
70
71 LOGICAL :: PARTICLE_SLIDE
72
73 LOGICAL, PARAMETER :: report_excess_overlap = .FALSE.
74
75
76
77
78 IF(USE_COHESION) PostCohesive(:) = ZERO
79
80 CALL CALC_DEM_FORCE_WITH_WALL_STL
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
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
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
132
133 (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
138 (LL) = PostCohesive(LL) + FORCE_COH
139 ENDIF
140 ENDIF
141 ENDIF
142
143 IF (ENERGY_EQ) THEN
144
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
149 (LL) = Q_Source(LL) + QQ_TMP
150
151
152 (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
172 = R_LM-DIST_MAG
173
174 IF (report_excess_overlap) call print_excess_overlap
175
176
177
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
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
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
202 (:) = -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
212 (:) = -KT_DES * PFT_TMP(:)
213 FTS2(:) = -ETAT_DES * V_REL_TANG
214 FT(:) = FTS1(:) + FTS2(:)
215
216
217
218 = .FALSE.
219 CALL CFSLIDE(V_REL_TANG(:), PARTICLE_SLIDE, MEW, &
220 FT(:), FN(:))
221
222
223
224 = 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
234
235 (:) = FC_TMP(:) + FN(:) + FT(:)
236
237 (1,LL) = FC(1,LL) + FC_TMP(1)
238
239 (2,LL) = FC(2,LL) + FC_TMP(2)
240
241 (3,LL) = FC(3,LL) + FC_TMP(3)
242
243 I = PAIRS(2,CC)
244
245 (1,I) = FC(1,I) - FC_TMP(1)
246
247 (2,I) = FC(2,I) - FC_TMP(2)
248
249 (3,I) = FC(3,I) - FC_TMP(3)
250
251
252
253
254 (1,LL) = TOW(1,LL) + TOW_TMP(1,1)
255
256 (2,LL) = TOW(2,LL) + TOW_TMP(2,1)
257
258 (3,LL) = TOW(3,LL) + TOW_TMP(3,1)
259
260
261 (1,I) = TOW(1,I) + TOW_TMP(1,2)
262
263 (2,I) = TOW(2,I) + TOW_TMP(2,2)
264
265 (3,I) = TOW(3,I) + TOW_TMP(3,2)
266
267
268
269 IF (PARTICLE_SLIDE) THEN
270
271
272 (:,CC) = -( FT(:) - FTS2(:) ) / KT_DES
273 ELSE
274 PFT_PAIR(:,CC) = PFT_TMP(:)
275 ENDIF
276
277 ENDDO
278
279
280
281
282 RETURN
283
284 contains
285
286
287
288
289
290
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
313
314
315
316
317
318
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
328 DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: PFT
329
330
331
332
333 DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: SIGMAT_OLD
334 LOGICAL, INTENT(IN) :: ALREADY_COLLIDING
335
336
337
338 DOUBLE PRECISION :: OVERLAP_T(3)
339
340 DOUBLE PRECISION, DIMENSION(3) :: SIGMAT
341
342
343
344 DOUBLE PRECISION :: DTSOLID_TMP
345
346
347
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
353
354 DOUBLE PRECISION :: TMP_AX(3), TMP_MAG
355
356
357
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
369
370
371
372
373
374
375 = 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
380 = des_crossprdct(tmp_ax,norm_old)
381
382 = 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
391 (:) = SIGMAT(:)
392 ELSE
393
394 = 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