File: /nfs/home/0/users/jenkins/mfix.git/model/calc_d_n.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 SUBROUTINE CALC_D_N(A_M, VXF_GS, VXF_SS, D_N, IER)
16
17
18
19
20 use physprop, only: MMAX
21
22 use run, only: MOMENTUM_Y_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(OUT) :: D_N(DIMENSION_3, 0:DIMENSION_M)
52
53 INTEGER, INTENT(OUT) :: IER
54
55
56
57
58
59
60 DOUBLE PRECISION, dimension(:,:), allocatable :: AM0
61
62 LOGICAL :: ANY_SOLIDS_Y_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_Y_MOMENTUM = any(MOMENTUM_Y_EQ(1:MMAX))
82
83 IF(MOMENTUM_Y_EQ(0))THEN
84 IF(ANY_SOLIDS_Y_MOMENTUM) THEN
85 CALL CALC_D_N_GAS_AND_SOLIDS(AM0, VxF_GS, VxF_SS, D_N)
86 ELSE
87 CALL CALC_D_N_GAS_ONLY(AM0, VxF_GS, D_N)
88 ENDIF
89
90 ELSEIF(ANY_SOLIDS_Y_MOMENTUM) THEN
91 CALL CALC_D_N_SOLIDS_ONLY(AM0, VxF_GS, VxF_SS, D_N)
92 ENDIF
93
94 deallocate(AM0)
95
96 RETURN
97 END SUBROUTINE CALC_D_N
98
99
100
101
102
103
104
105
106
107
108
109
110 SUBROUTINE CALC_D_N_GAS_AND_SOLIDS(AM0, VXF_GS, VXF_SS, D_N)
111
112
113
114
115 use fldvar, only: EP_G, EP_s
116
117 use physprop, only: MMAX
118
119 use run, only: MOMENTUM_Y_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: AXZ
128
129 use fun_avg, only: AVG_Y
130
131 use compar, only: IJKStart3, IJKEnd3
132
133 use functions, only: FUNLM, NORTH_OF
134
135 use functions, only: IP_AT_N, MFLOW_AT_N
136
137 use indices, only: J_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_N(DIMENSION_3, 0:DIMENSION_M)
161
162
163
164
165 INTEGER :: LM, M, L, LPL, LP
166 INTEGER :: J, IJK, IJKN
167
168 DOUBLE PRECISION :: EPGA
169
170 DOUBLE PRECISION :: EPSA(DIMENSION_M)
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 IF (IP_AT_N(IJK) .OR. MFLOW_AT_N(IJK)) THEN
209 D_N(IJK,0:MMAX) = ZERO
210 CYCLE
211 ENDIF
212
213 AREA_FACE = merge(ONE, AXZ(IJK), CARTESIAN_GRID)
214 J = J_OF(IJK)
215 IJKN = NORTH_OF(IJK)
216 EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
217
218 SUM_VXF_GS = ZERO
219 DO M= 1, MMAX
220 EPSA(M) = AVG_Y(EP_S(IJK,M),EP_S(IJKN,M),J)
221
222 = SUM_VXF_GS + VXF_GS(IJK,M)
223 SUM_VXF_SS(M) = ZERO
224 DO L = 1, MMAX
225 IF (L .NE. M) THEN
226 LM = FUNLM(L,M)
227
228 (M) = SUM_VXF_SS(M) + VXF_SS(IJK,LM)
229 ENDIF
230 ENDDO
231 ENDDO
232
233
234 = ZERO
235 DEN_MGas = ZERO
236 DO M= 1, MMAX
237 IF(MOMENTUM_Y_EQ(M)) THEN
238 NUM_MGas = NUM_MGas + (EPSA(M)*VXF_GS(IJK,M) / &
239 ((-AM0(IJK,M)) + VXF_GS(IJK,M) + SUM_VXF_SS(M) + &
240 SMALL_NUMBER))
241 DEN_MGas = DEN_MGas + (VXF_GS(IJK,M)* &
242 ((-AM0(IJK,M)) + SUM_VXF_SS(M)) / &
243 ((-AM0(IJK,M)) + VXF_GS(IJK,M) + SUM_VXF_SS(M) + &
244 SMALL_NUMBER))
245 ELSE
246 DEN_MGas = DEN_MGas + VXF_GS(IJK,M)
247 ENDIF
248 ENDDO
249
250
251
252 IF (MODEL_B) THEN
253
254 = -AM0(IJK,0) + DEN_MGas
255 IF(abs(TMPdp) > SMALL_NUMBER ) THEN
256 D_N(IJK,0) = P_SCALE*AREA_FACE/TMPdp
257 ELSE
258 D_N(IJK,0) = ZERO
259 ENDIF
260
261 DO M = 1, MMAX
262 IF(MOMENTUM_Y_EQ(M)) THEN
263 TMPdp = -AM0(IJK,M) + VXF_GS(IJK,M)
264 IF(abs(TMPdp) > SMALL_NUMBER) THEN
265 D_N(IJK,M) = D_N(IJK,0)*VXF_GS(IJK,M) / TMPdp
266 ELSE
267 D_N(IJK,M) = ZERO
268 ENDIF
269 ELSE
270 D_N(IJK,M) = ZERO
271 ENDIF
272 ENDDO
273
274
275
276 ELSE
277
278
279 = -AM0(IJK,0) + DEN_MGas
280 IF(abs(TMPdp) > SMALL_NUMBER) THEN
281 D_N(IJK,0) = P_SCALE*AREA_FACE*(EPGA+NUM_MGas) / TMPdp
282 ELSE
283 D_N(IJK,0) = ZERO
284 ENDIF
285
286 DO M= 1, MMAX
287
288 = VXF_GS(IJK,M)*EPGA/ &
289 ((-AM0(IJK,0))+SUM_VXF_GS+SMALL_NUMBER)
290 DEN_MSol_LGas = VXF_GS(IJK,M)*( &
291 ((-AM0(IJK,0))+(SUM_VXF_GS - VXF_GS(IJK,M))) / &
292 ((-AM0(IJK,0))+SUM_VXF_GS + SMALL_NUMBER) )
293
294
295 (M) = ZERO
296 DEN_MSol_LSol(M) = ZERO
297 DO L = 1, MMAX
298 IF (L .NE. M) THEN
299 LM = FUNLM(L,M)
300 SUM_VXF_SS_wt_M = ZERO
301 DO Lp = 1, MMAX
302 IF((Lp /= L) .AND. (Lp /= M)) THEN
303 LpL = FUNLM(Lp,L)
304
305 = SUM_VXF_SS_wt_M + &
306 VXF_SS(IJK,LpL)
307 ENDIF
308 ENDDO
309
310 IF(MOMENTUM_Y_EQ(L)) THEN
311 NUM_MSol_LSol(M) = NUM_MSol_LSol(M) + &
312 (VXF_SS(IJK,LM)*EPSA(L) / &
313 ((-AM0(IJK,L))+VXF_GS(IJK,L)+SUM_VXF_SS(L)+ &
314 SMALL_NUMBER ))
315
316 DEN_MSol_LSol(M) = DEN_MSol_LSol(M) + &
317 VXF_SS(IJK,LM)*(((-AM0(IJK,L)) + &
318 VXF_GS(IJK,L) + SUM_VXF_SS_wt_M ) / &
319 ((-AM0(IJK,L)) + VXF_GS(IJK,L)+ &
320 SUM_VXF_SS(L) + SMALL_NUMBER) )
321 ELSE
322 DEN_MSol_LSol(M) = DEN_MSol_LSol(M) + &
323 VXF_SS(IJK,LM)
324 ENDIF
325 ENDIF
326 ENDDO
327
328
329 IF(MOMENTUM_Y_EQ(M)) THEN
330 TMPdp = -AM0(IJK,M) + DEN_MSol_LGas + DEN_MSol_LSol(M)
331
332 IF(abs(TMPdp) > SMALL_NUMBER) THEN
333 D_N(IJK,M) = P_SCALE*AREA_FACE*(EPSA(M) + &
334 NUM_MSol_LSol(M) + NUM_MSol_LGas) / TMPdp
335 ELSE
336 D_N(IJK,M) = ZERO
337 ENDIF
338 ELSE
339 D_N(IJK,M) = ZERO
340 ENDIF
341 ENDDO
342
343 ENDIF
344 ENDDO
345
346
347 RETURN
348 END SUBROUTINE CALC_D_N_GAS_AND_SOLIDS
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364 SUBROUTINE CALC_D_N_GAS_ONLY(AM0, VXF_GS, D_N)
365
366
367
368
369 use fldvar, only: EP_G
370
371 use physprop, only: MMAX
372
373 use run, only: MODEL_B
374
375 use scales, only: P_SCALE
376
377 use cutcell, only: CARTESIAN_GRID
378
379 use geometry, only: AXZ
380
381 use geometry, only: VOL_V
382
383 use fun_avg, only: AVG_Y
384
385 use compar, only: IJKStart3, IJKEnd3
386
387 use functions, only: IP_AT_N, MFLOW_AT_N, NORTH_OF
388
389 use indices, only: J_OF
390
391 use qmom_kinetic_equation, only: QMOMK, QMOMK_NN
392 use qmom_kinetic_equation, only: QMOMK_F_GS, QMOMK_F_GS
393
394
395
396
397 use param, only: DIMENSION_3, DIMENSION_M
398
399 use param1, only: DIMENSION_LM
400
401 use param1, only: ZERO, SMALL_NUMBER, ONE
402
403 IMPLICIT NONE
404
405
406
407
408 DOUBLE PRECISION, INTENT(IN) :: AM0(DIMENSION_3, 0:DIMENSION_M)
409
410 DOUBLE PRECISION, INTENT(IN) :: VxF_gs(DIMENSION_3, DIMENSION_M)
411
412
413 DOUBLE PRECISION, INTENT(INOUT) :: d_n(DIMENSION_3, 0:DIMENSION_M)
414
415
416
417
418 INTEGER :: M
419 INTEGER :: INN
420 INTEGER :: J, IJK, IJKN
421
422 DOUBLE PRECISION :: EPGA
423
424 DOUBLE PRECISION :: AREA_FACE
425
426 DOUBLE PRECISION :: SUM_VXF_GS
427
428 DOUBLE PRECISION :: TMPdp
429
430
431
432
433
434
435
436 DO IJK = IJKSTART3, IJKEND3
437
438 IF (IP_AT_N(IJK) .OR. MFLOW_AT_N(IJK)) THEN
439 D_N(IJK,0) = ZERO
440 CYCLE
441 ENDIF
442
443 AREA_FACE = merge(ONE, AXZ(IJK), CARTESIAN_GRID)
444 J = J_OF(IJK)
445 IJKN = NORTH_OF(IJK)
446 EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
447
448 SUM_VXF_GS = ZERO
449 IF (.NOT. QMOMK) THEN
450 DO M= 1, MMAX
451
452 = 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_V(IJK)*AVG_Y( &
458 QMOMK_F_GS(INN,IJK,M),QMOMK_F_GS(INN,IJKN,M),J)
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_N(IJK,0) = P_SCALE*AREA_FACE/TMPdp
467 ELSE
468 D_N(IJK,0) = P_SCALE*AREA_FACE*EPGA/TMPdp
469 ENDIF
470 ELSE
471 D_N(IJK,0) = ZERO
472 ENDIF
473
474 ENDDO
475
476
477 RETURN
478 END SUBROUTINE CALC_D_N_GAS_ONLY
479
480
481
482
483
484
485
486
487
488
489
490
491
492 SUBROUTINE CALC_D_N_SOLIDS_ONLY(AM0, VXF_GS, VXF_SS, D_N)
493
494
495
496
497 use fldvar, only: EP_G, EP_s
498
499 use physprop, only: MMAX
500
501 use run, only: MOMENTUM_Y_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: AXZ
510
511 use fun_avg, only: AVG_Y
512
513 use compar, only: IJKStart3, IJKEnd3
514
515 use functions, only: FUNLM, NORTH_OF
516
517 use functions, only: IP_AT_N, MFLOW_AT_N
518
519 use indices, only: J_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_N(DIMENSION_3, 0:DIMENSION_M)
543
544
545
546
547 INTEGER :: LM, M, L, LPL, LP
548 INTEGER :: J, IJK, IJKN
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_N(IJK) .OR. MFLOW_AT_N(IJK) .OR. MODEL_B) THEN
576 D_N(IJK,1:MMAX) = ZERO
577 CYCLE
578 ENDIF
579
580 AREA_FACE = merge(ONE, AXZ(IJK), CARTESIAN_GRID)
581 J = J_OF(IJK)
582 IJKN = NORTH_OF(IJK)
583
584 DO M= 1, MMAX
585 EPSA(M) = AVG_Y(EP_S(IJK,M),EP_S(IJKN,M),J)
586 SUM_VXF_SS(M) = ZERO
587 DO L = 1, MMAX
588 IF (L .NE. M) THEN
589 LM = FUNLM(L,M)
590
591 (M) = SUM_VXF_SS(M) + VXF_SS(IJK,LM)
592 ENDIF
593 ENDDO
594 ENDDO
595
596 DO M= 1, MMAX
597
598 (M) = ZERO
599 DEN_MSol_LSol(M) = ZERO
600 DO L = 1, MMAX
601 IF (L .NE. M) THEN
602 LM = FUNLM(L,M)
603 SUM_VXF_SS_wt_M = ZERO
604 DO Lp = 1, MMAX
605 IF((Lp /= L) .AND. (Lp /= M)) THEN
606 LpL = FUNLM(Lp,L)
607
608 =SUM_VXF_SS_wt_M+VXF_SS(IJK,LpL)
609 ENDIF
610 ENDDO
611 IF(MOMENTUM_Y_EQ(L)) THEN
612 NUM_MSol_LSol(M) = NUM_MSol_LSol(M) + &
613 (VXF_SS(IJK,LM)*EPSA(L)/((-AM0(IJK,L)) + &
614 VXF_GS(IJK,L)+SUM_VXF_SS(L)+SMALL_NUMBER ))
615 DEN_MSol_LSol(M) = DEN_MSol_LSol(M) + &
616 VXF_SS(IJK,LM)*(((-AM0(IJK,L)) + VXF_GS(IJK,L)+&
617 SUM_VXF_SS_wt_M )/((-AM0(IJK,L))+VXF_GS(IJK,L)+&
618 SUM_VXF_SS(L) + SMALL_NUMBER))
619 ELSE
620 DEN_MSol_LSol(M) = DEN_MSol_LSol(M)+VXF_SS(IJK,LM)
621 ENDIF
622 ENDIF
623 ENDDO
624
625
626 IF (MOMENTUM_Y_EQ(M)) THEN
627 TMPdp = -AM0(IJK,M) + VXF_GS(IJK,M) + DEN_MSol_LSol(M)
628 IF(abs(TMPdp) > SMALL_NUMBER) THEN
629 D_N(IJK,M) = P_SCALE*AREA_FACE* (EPSA(M) + &
630 NUM_MSol_LSol(M))/TMPdp
631 ELSE
632 D_N(IJK,M) = ZERO
633 ENDIF
634 ELSE
635 D_N(IJK,M) = ZERO
636 ENDIF
637 ENDDO
638 ENDDO
639
640
641 RETURN
642 END SUBROUTINE CALC_D_N_SOLIDS_ONLY
643