File: /nfs/home/0/users/jenkins/mfix.git/model/calc_d_t.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 SUBROUTINE CALC_D_T(A_M, VXF_GS, VXF_SS, D_T, IER)
16
17
18
19
20 use physprop, only: MMAX
21
22 use run, only: MOMENTUM_Z_EQ
23
24 use discretelement, only: DES_CONTINUUM_COUPLED
25
26 use discretelement, only: DES_CONTINUUM_HYBRID
27
28 use discretelement, only: VXF_GDS, VXF_SDS
29
30 use physprop, only: MMAX
31
32
33
34
35 use param, only: DIMENSION_3, DIMENSION_M
36
37 use param1, only: DIMENSION_LM
38
39 IMPLICIT NONE
40
41
42
43
44 DOUBLE PRECISION, INTENT(IN):: A_m(DIMENSION_3,-3:3,0:DIMENSION_M)
45
46 DOUBLE PRECISION, INTENT(IN) :: VxF_gs(DIMENSION_3, DIMENSION_M)
47
48 DOUBLE PRECISION, INTENT(IN) :: VxF_ss(DIMENSION_3, DIMENSION_LM)
49
50
51 DOUBLE PRECISION, INTENT(INOUT) :: d_t(DIMENSION_3, 0:DIMENSION_M)
52
53 INTEGER, INTENT(INOUT) :: IER
54
55
56
57
58
59
60 DOUBLE PRECISION, dimension(:,:), allocatable :: AM0
61
62 LOGICAL :: ANY_SOLIDS_Z_MOMENTUM
63
64
65
66 allocate(AM0(DIMENSION_3, 0:DIMENSION_M))
67
68
69 = 0
70
71
72 (:,:) = A_M(:,0,:)
73
74
75 IF(DES_CONTINUUM_COUPLED) THEN
76 AM0(:,0) = AM0(:,0) - VXF_GDS(:)
77 IF (DES_CONTINUUM_HYBRID) &
78 AM0(:,1:MMAX) = AM0(:,1:MMAX) - VXF_SDS(:,1:MMAX)
79 ENDIF
80
81 ANY_SOLIDS_Z_MOMENTUM = any(MOMENTUM_Z_EQ(1:MMAX))
82
83
84 IF (MOMENTUM_Z_EQ(0)) THEN
85 IF(ANY_SOLIDS_Z_MOMENTUM)THEN
86 CALL CALC_D_T_GAS_AND_SOLIDS(AM0, VXF_GS, VXF_SS, D_T)
87 ELSE
88 CALL CALC_D_T_GAS_ONLY(AM0, VXF_GS, D_T)
89 ENDIF
90 ELSEIF (ANY_SOLIDS_Z_MOMENTUM) THEN
91 CALL CALC_D_T_SOLIDS_ONLY(AM0, VXF_GS, VXF_SS, D_T)
92 ENDIF
93
94 deallocate(AM0)
95
96 RETURN
97 END SUBROUTINE CALC_D_T
98
99
100
101
102
103
104
105
106
107
108
109
110 SUBROUTINE CALC_D_T_GAS_AND_SOLIDS(AM0, VXF_GS, VXF_SS, D_T)
111
112
113
114
115 use fldvar, only: EP_G, EP_s
116
117 use physprop, only: MMAX
118
119 use run, only: MOMENTUM_Z_EQ
120
121 use run, only: MODEL_B
122
123 use scales, only: P_SCALE
124
125 use cutcell, only: CARTESIAN_GRID
126
127 use geometry, only: AXY
128
129 use fun_avg, only: AVG_Z
130
131 use compar, only: IJKStart3, IJKEnd3
132
133 use functions, only: FUNLM, TOP_OF
134
135 use functions, only: IP_AT_T, MFLOW_AT_T
136
137 use indices, only: K_OF
138
139
140
141
142 use param, only: DIMENSION_3, DIMENSION_M
143
144 use param1, only: DIMENSION_LM
145
146 use param1, only: ZERO, SMALL_NUMBER, ONE
147
148 IMPLICIT NONE
149
150
151
152
153 DOUBLE PRECISION, INTENT(IN) :: AM0(DIMENSION_3, 0:DIMENSION_M)
154
155 DOUBLE PRECISION, INTENT(IN) :: VxF_gs(DIMENSION_3, DIMENSION_M)
156
157 DOUBLE PRECISION, INTENT(IN) :: VxF_ss(DIMENSION_3, DIMENSION_LM)
158
159
160 DOUBLE PRECISION, INTENT(INOUT) :: D_T(DIMENSION_3, 0:DIMENSION_M)
161
162
163
164
165 INTEGER :: LM, M, L, LPL, LP
166 INTEGER :: IJK, IJKT, K
167
168 DOUBLE PRECISION :: EPSA(DIMENSION_M)
169
170 DOUBLE PRECISION :: EPGA
171
172 DOUBLE PRECISION :: AREA_FACE
173
174 DOUBLE PRECISION :: SUM_VXF_GS
175
176 DOUBLE PRECISION :: SUM_VXF_SS(DIMENSION_M)
177
178 DOUBLE PRECISION :: SUM_VXF_SS_wt_M
179
180
181 DOUBLE PRECISION :: DEN_MGas, NUM_MGas
182
183
184
185 DOUBLE PRECISION :: NUM_MSol_LSol(DIMENSION_M)
186
187
188
189 DOUBLE PRECISION :: DEN_MSol_LSol(DIMENSION_M)
190
191
192
193 DOUBLE PRECISION :: NUM_MSol_LGas
194
195
196
197 DOUBLE PRECISION :: DEN_MSol_LGas
198
199 DOUBLE PRECISION :: TMPdp
200
201
202
203
204
205
206
207 DO IJK = IJKSTART3, IJKEND3
208
209 IF (IP_AT_T(IJK) .OR. MFLOW_AT_T(IJK)) THEN
210 DO M= 0, MMAX
211 D_T(IJK,M) = ZERO
212 ENDDO
213 CYCLE
214 ENDIF
215
216 AREA_FACE = merge(ONE, AXY(IJK), CARTESIAN_GRID)
217 K = K_OF(IJK)
218 IJKT = TOP_OF(IJK)
219 EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
220
221 SUM_VXF_GS = ZERO
222 DO M= 1, MMAX
223 EPSA(M) = AVG_Z(EP_S(IJK,M),EP_S(IJKT,M),K)
224
225 = SUM_VXF_GS + VXF_GS(IJK,M)
226 SUM_VXF_SS(M) = ZERO
227 DO L = 1, MMAX
228 IF (L .NE. M) THEN
229 LM = FUNLM(L,M)
230
231 (M) = SUM_VXF_SS(M) + VXF_SS(IJK,LM)
232 ENDIF
233 ENDDO
234 ENDDO
235
236
237 = ZERO
238 NUM_MGas = ZERO
239 DO M= 1, MMAX
240 IF(MOMENTUM_Z_EQ(M)) THEN
241 NUM_MGas = NUM_MGas + (EPSA(M)*VXF_GS(IJK,M)/ &
242 ((-AM0(IJK,M))+VXF_GS(IJK,M)+SUM_VXF_SS(M) + &
243 SMALL_NUMBER))
244 DEN_MGas = DEN_MGas + (VXF_GS(IJK,M)* &
245 ((-AM0(IJK,M))+SUM_VXF_SS(M))/ &
246 ((-AM0(IJK,M))+VXF_GS(IJK,M)+SUM_VXF_SS(M)+ &
247 SMALL_NUMBER))
248 ELSE
249 DEN_MGas = DEN_MGas + VXF_GS(IJK,M)
250 ENDIF
251 ENDDO
252
253
254
255 IF (MODEL_B) THEN
256
257 = -AM0(IJK,0) + DEN_MGas
258 IF(abs(TMPdp) > SMALL_NUMBER) THEN
259 D_T(IJK,0) = P_SCALE*AREA_FACE/TMPdp
260 ELSE
261 D_T(IJK,0) = ZERO
262 ENDIF
263
264 DO M = 1, MMAX
265 IF(MOMENTUM_Z_EQ(M)) THEN
266 TMPdp = -AM0(IJK,M) + VXF_GS(IJK,M)
267 IF(abs(TMPdp) > SMALL_NUMBER) THEN
268 D_T(IJK,M) = D_T(IJK,0)*VXF_GS(IJK,M)/TMPdp
269 ELSE
270 D_T(IJK,M) = ZERO
271 ENDIF
272 ELSE
273 D_T(IJK,M) = ZERO
274 ENDIF
275 ENDDO
276
277
278
279 ELSE
280
281
282 = -AM0(IJK,0)+DEN_MGas
283 IF(abs(TMPdp) > SMALL_NUMBER) THEN
284 D_T(IJK,0) = P_SCALE*AREA_FACE*(EPGA+NUM_MGas)/TMPdp
285 ELSE
286 D_T(IJK,0) = ZERO
287 ENDIF
288
289 DO M= 1, MMAX
290
291 = VXF_GS(IJK,M)*EPGA/ &
292 ((-AM0(IJK,0))+SUM_VXF_GS+SMALL_NUMBER)
293 DEN_MSol_LGas = VXF_GS(IJK,M)*( &
294 ((-AM0(IJK,0)) + (SUM_VXF_GS - VXF_GS(IJK,M)))/ &
295 ((-AM0(IJK,0)) + SUM_VXF_GS+SMALL_NUMBER))
296
297 (M) = ZERO
298 DEN_MSol_LSol(M) = ZERO
299 DO L = 1, MMAX
300 IF (L .NE. M) THEN
301 LM = FUNLM(L,M)
302 SUM_VXF_SS_wt_M = ZERO
303 DO Lp = 1, MMAX
304 IF( (Lp /= L) .AND. (Lp /= M) ) THEN
305 LpL = FUNLM(Lp,L)
306
307 = &
308 SUM_VXF_SS_wt_M + VXF_SS(IJK,LpL)
309 ENDIF
310 ENDDO
311
312 IF(MOMENTUM_Z_EQ(L)) THEN
313 NUM_MSol_LSol(M) = NUM_MSol_LSol(M) + &
314 (VXF_SS(IJK,LM)*EPSA(L)/((-AM0(IJK,L)) + &
315 VXF_GS(IJK,L)+SUM_VXF_SS(L)+ SMALL_NUMBER))
316
317 DEN_MSol_LSol(M) = DEN_MSol_LSol(M) + &
318 VXF_SS(IJK,LM)*(((-AM0(IJK,L)) + &
319 VXF_GS(IJK,L)+SUM_VXF_SS_wt_M )/ &
320 ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS(L)+ &
321 SMALL_NUMBER))
322 ELSE
323 DEN_MSol_LSol(M) = &
324 DEN_MSol_LSol(M) + VXF_SS(IJK,LM)
325 ENDIF
326 ENDIF
327 ENDDO
328
329
330 IF(MOMENTUM_Z_EQ(M)) THEN
331 TMPdp = -AM0(IJK,M) + DEN_MSol_LGas + DEN_MSol_LSol(M)
332 IF(abs(TMPdp) > SMALL_NUMBER) THEN
333 D_T(IJK,M) = P_SCALE*AREA_FACE*(EPSA(M) + &
334 NUM_MSol_LSol(M) + NUM_MSol_LGas)/TMPdp
335 ELSE
336 D_T(IJK,M) = ZERO
337 ENDIF
338 ELSE
339 D_T(IJK,M) = ZERO
340 ENDIF
341 ENDDO
342
343 ENDIF
344 ENDDO
345
346
347 RETURN
348 END SUBROUTINE CALC_D_T_GAS_AND_SOLIDS
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365 SUBROUTINE CALC_D_T_GAS_ONLY(AM0, VXF_GS, D_T)
366
367
368
369
370 use fldvar, only: EP_G
371
372 use physprop, only: MMAX
373
374 use run, only: MODEL_B
375
376 use scales, only: P_SCALE
377
378 use cutcell, only: CARTESIAN_GRID
379
380 use geometry, only: AXY
381
382 use geometry, only: VOL_W
383
384 use fun_avg, only: AVG_Z
385
386 use compar, only: IJKStart3, IJKEnd3
387
388 use functions, only: IP_AT_T, MFLOW_AT_T, TOP_OF
389
390 use indices, only: K_OF
391
392 use qmom_kinetic_equation, only: QMOMK, QMOMK_NN
393 use qmom_kinetic_equation, only: QMOMK_F_GS, QMOMK_F_GS
394
395
396
397
398 use param, only: DIMENSION_3, DIMENSION_M
399
400 use param1, only: DIMENSION_LM
401
402 use param1, only: ZERO, SMALL_NUMBER, ONE
403
404 IMPLICIT NONE
405
406
407
408
409 DOUBLE PRECISION, INTENT(IN) :: AM0(DIMENSION_3, 0:DIMENSION_M)
410
411 DOUBLE PRECISION, INTENT(IN) :: VxF_gs(DIMENSION_3, DIMENSION_M)
412
413
414 DOUBLE PRECISION, INTENT(INOUT) :: d_t(DIMENSION_3, 0:DIMENSION_M)
415
416
417
418
419
420 INTEGER :: M
421 INTEGER :: INN
422 INTEGER :: IJK, IJKT, K
423
424 DOUBLE PRECISION :: EPGA
425
426 DOUBLE PRECISION :: AREA_FACE
427
428 DOUBLE PRECISION :: SUM_VXF_GS
429
430 DOUBLE PRECISION :: TMPdp
431
432
433
434
435
436
437
438 DO IJK = ijkstart3, ijkend3
439 IF (IP_AT_T(IJK) .OR. MFLOW_AT_T(IJK)) THEN
440 D_T(IJK,0) = ZERO
441 CYCLE
442 ENDIF
443
444 AREA_FACE = merge(ONE, AXY(IJK), CARTESIAN_GRID)
445 K = K_OF(IJK)
446 IJKT = TOP_OF(IJK)
447 EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
448
449 SUM_VXF_GS = ZERO
450 IF (.NOT. QMOMK) THEN
451 DO M= 1, MMAX
452 SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)
453 ENDDO
454 ELSE
455 DO INN = 1, QMOMK_NN
456 DO M = 1, MMAX
457 SUM_VXF_GS = SUM_VXF_GS + VOL_W(IJK)*AVG_Z( &
458 QMOMK_F_GS(INN,IJK,M), QMOMK_F_GS(INN,IJKT,M),K)
459 ENDDO
460 ENDDO
461 ENDIF
462
463 TMPdp = -AM0(IJK,0)+SUM_VXF_GS
464 IF(abs(TMPdp) > SMALL_NUMBER) THEN
465 IF (MODEL_B) THEN
466 D_T(IJK,0) = P_SCALE*AREA_FACE/TMPdp
467 ELSE
468 D_T(IJK,0) = P_SCALE*AREA_FACE*EPGA/TMPdp
469 ENDIF
470 ELSE
471 D_T(IJK,0) = ZERO
472 ENDIF
473 ENDDO
474
475
476 RETURN
477 END SUBROUTINE CALC_D_T_GAS_ONLY
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492 SUBROUTINE CALC_D_T_SOLIDS_ONLY(AM0, VXF_GS, VXF_SS, D_T)
493
494
495
496
497 use fldvar, only: EP_G, EP_s
498
499 use physprop, only: MMAX
500
501 use run, only: MOMENTUM_Z_EQ
502
503 use run, only: MODEL_B
504
505 use scales, only: P_SCALE
506
507 use cutcell, only: CARTESIAN_GRID
508
509 use geometry, only: AXY
510
511 use fun_avg, only: AVG_Z
512
513 use compar, only: IJKStart3, IJKEnd3
514
515 use functions, only: FUNLM, TOP_OF
516
517 use functions, only: IP_AT_T, MFLOW_AT_T
518
519 use indices, only: K_OF
520
521
522
523
524 use param, only: DIMENSION_3, DIMENSION_M
525
526 use param1, only: DIMENSION_LM
527
528 use param1, only: ZERO, SMALL_NUMBER, ONE
529
530 IMPLICIT NONE
531
532
533
534
535 DOUBLE PRECISION, INTENT(IN) :: AM0(DIMENSION_3, 0:DIMENSION_M)
536
537 DOUBLE PRECISION, INTENT(IN) :: VxF_gs(DIMENSION_3, DIMENSION_M)
538
539 DOUBLE PRECISION, INTENT(IN) :: VxF_ss(DIMENSION_3, DIMENSION_LM)
540
541
542 DOUBLE PRECISION, INTENT(INOUT) :: d_t(DIMENSION_3, 0:DIMENSION_M)
543
544
545
546
547 INTEGER :: LM, M, L, LPL, LP
548 INTEGER :: IJK, IJKT, K
549
550 DOUBLE PRECISION :: EPSA(DIMENSION_M)
551
552 DOUBLE PRECISION :: AREA_FACE
553
554 DOUBLE PRECISION :: SUM_VXF_SS(DIMENSION_M)
555
556 DOUBLE PRECISION :: SUM_VXF_SS_wt_M
557
558
559
560 DOUBLE PRECISION :: NUM_MSol_LSol(DIMENSION_M)
561
562
563
564 DOUBLE PRECISION :: DEN_MSol_LSol(DIMENSION_M)
565
566 DOUBLE PRECISION :: TMPdp
567
568
569
570
571
572
573
574 DO IJK = ijkstart3, ijkend3
575 IF (IP_AT_T(IJK) .OR. MFLOW_AT_T(IJK) .OR. MODEL_B) THEN
576 DO M= 1, MMAX
577 D_T(IJK,M) = ZERO
578 ENDDO
579 CYCLE
580 ENDIF
581
582 AREA_FACE = merge(ONE, AXY(IJK), CARTESIAN_GRID)
583 K = K_OF(IJK)
584 IJKT = TOP_OF(IJK)
585
586 DO M= 1, MMAX
587 EPSA(M) = AVG_Z(EP_S(IJK,M),EP_S(IJKT,M),K)
588 SUM_VXF_SS(M) = ZERO
589 DO L = 1, MMAX
590 IF (L .NE. M) THEN
591 LM = FUNLM(L,M)
592
593 (M) = SUM_VXF_SS(M) + VXF_SS(IJK,LM)
594 ENDIF
595 ENDDO
596 ENDDO
597
598 DO M= 1, MMAX
599
600 (M) = ZERO
601 DEN_MSol_LSol(M) = ZERO
602 DO L = 1, MMAX
603 IF (L .NE. M) THEN
604 LM = FUNLM(L,M)
605 SUM_VXF_SS_wt_M = ZERO
606 DO Lp = 1, MMAX
607 IF ( (Lp .NE. L) .AND. (Lp .NE. M) ) THEN
608 LpL = FUNLM(Lp,L)
609
610 = SUM_VXF_SS_wt_M + VXF_SS(IJK,LpL)
611 ENDIF
612 ENDDO
613 IF (MOMENTUM_Z_EQ(L)) THEN
614 NUM_MSol_LSol(M) = NUM_MSol_LSol(M) + &
615 (VXF_SS(IJK,LM)*EPSA(L)/&
616 ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS(L)+&
617 SMALL_NUMBER))
618 DEN_MSol_LSol(M) = DEN_MSol_LSol(M) +VXF_SS(IJK,LM)*(&
619 ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS_wt_M )/&
620 ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS(L)+&
621 SMALL_NUMBER ))
622 ELSE
623 DEN_MSol_LSol(M) = DEN_MSol_LSol(M) + VXF_SS(IJK,LM)
624 ENDIF
625 ENDIF
626 ENDDO
627
628
629 IF(MOMENTUM_Z_EQ(M)) THEN
630 TMPdp = (-AM0(IJK,M))+VXF_GS(IJK,M)+DEN_MSol_LSol(M)
631
632 IF(abs(TMPdp) > SMALL_NUMBER) THEN
633 D_T(IJK,M) = P_SCALE*AREA_FACE*&
634 (EPSA(M) + NUM_MSol_LSol(M))/TMPdp
635 ELSE
636 D_T(IJK,M) = ZERO
637 ENDIF
638 ELSE
639 D_T(IJK,M) = ZERO
640 ENDIF
641 ENDDO
642 ENDDO
643
644
645 RETURN
646 END SUBROUTINE CALC_D_T_SOLIDS_ONLY
647