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