File: /nfs/home/0/users/jenkins/mfix.git/model/set_outflow.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35 SUBROUTINE SET_OUTFLOW(BCV)
36
37
38
39 use bc, only: bc_type, bc_plane
40 use bc, only: bc_k_b, bc_k_t
41 use bc, only: bc_j_s, bc_j_n
42 use bc, only: bc_i_w, bc_i_e
43
44 use run, only: kt_type_enum, ghd_2007
45 use fldvar, only: rop_g, rop_s
46 use fldvar, only: u_g, v_g, w_g
47 use fldvar, only: u_s, v_s, w_s
48 use physprop, only: mmax
49
50 use functions, only: im_of, ip_of, jm_of, jp_of, km_of, kp_of
51 use functions, only: fluid_at
52 use functions, only: is_on_mype_plus2layers
53 use functions, only: funijk
54 use compar, only: dead_cell_at
55
56 use param, only: dimension_m
57 use param1, only: undefined, zero
58 IMPLICIT NONE
59
60
61
62
63 INTEGER, INTENT(IN) :: BCV
64
65
66
67
68 INTEGER :: I, J, K, M
69
70 INTEGER :: IJK
71
72 INTEGER :: FIJK
73
74
75 DOUBLE PRECISION :: RVEL_G, RVEL_S(DIMENSION_M)
76
77
78
79 DO K = BC_K_B(BCV), BC_K_T(BCV)
80 DO J = BC_J_S(BCV), BC_J_N(BCV)
81 DO I = BC_I_W(BCV), BC_I_E(BCV)
82
83 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
84 IF(DEAD_CELL_AT(I,J,K)) CYCLE
85 IJK = FUNIJK(I,J,K)
86
87
88
89 IF (FLUID_AT(IM_OF(IJK))) THEN
90 FIJK = IM_OF(IJK)
91 RVEL_G = U_G(FIJK)
92 IF (MMAX>0) RVEL_S(:MMAX) = U_S(FIJK,:MMAX)
93
94 CALL SET_OUTFLOW_MISC(BCV, IJK, FIJK)
95 CALL SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
96 IF (BC_TYPE(BCV)=='P_INFLOW') &
97 CALL SET_PINOUTFLOW(BCV, IJK, FIJK, RVEL_G, RVEL_S)
98
99
100
101
102
103
104
105
106
107
108
109 IF (ROP_G(IJK) > ZERO) THEN
110 U_G(IJK) = ROP_G(FIJK)*U_G(FIJK)/ROP_G(IJK)
111 ELSE
112 U_G(IJK) = ZERO
113 ENDIF
114
115
116
117 (IJK) = V_G(FIJK)
118 W_G(IJK) = W_G(FIJK)
119
120 IF (MMAX > 0) THEN
121 WHERE (ROP_S(IJK,:MMAX) > ZERO)
122 U_S(IJK,:MMAX) = ROP_S(FIJK,:MMAX)*&
123 U_S(FIJK,:MMAX)/ROP_S(IJK,:MMAX)
124 ELSEWHERE
125 U_S(IJK,:MMAX) = ZERO
126 END WHERE
127 V_S(IJK,:MMAX) = V_S(FIJK,:MMAX)
128 W_S(IJK,:MMAX) = W_S(FIJK,:MMAX)
129 ENDIF
130
131 CALL SET_OUTFLOW_FLUXES(IJK, FIJK)
132 ENDIF
133
134
135
136
137 IF (FLUID_AT(IP_OF(IJK))) THEN
138 FIJK = IP_OF(IJK)
139
140
141 = -U_G(IJK)
142 IF (MMAX >0) RVEL_S(:MMAX) = -U_S(IJK,:MMAX)
143
144 CALL SET_OUTFLOW_MISC(BCV, IJK, FIJK)
145 CALL SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
146 IF (BC_TYPE(BCV)=='P_INFLOW') &
147 CALL SET_PINOUTFLOW(BCV, IJK, FIJK, RVEL_G, RVEL_S)
148
149
150
151
152
153
154 IF (U_G(IJK) == UNDEFINED) THEN
155 IF (ROP_G(IJK) > ZERO) THEN
156 U_G(IJK) = ROP_G(FIJK)*U_G(FIJK)/ROP_G(IJK)
157 ELSE
158 U_G(IJK) = ZERO
159 ENDIF
160 ENDIF
161 V_G(IJK) = V_G(FIJK)
162 W_G(IJK) = W_G(FIJK)
163
164 DO M = 1, MMAX
165 IF (U_S(IJK,M) == UNDEFINED) THEN
166 IF (ROP_S(IJK,M) > ZERO) THEN
167 U_S(IJK,M) = ROP_S(FIJK,M)*&
168 U_S(FIJK,M)/ROP_S(IJK,M)
169 ELSE
170 U_S(IJK,M) = ZERO
171 ENDIF
172 ENDIF
173 V_S(IJK,M) = V_S(FIJK,M)
174 W_S(IJK,M) = W_S(FIJK,M)
175 ENDDO
176
177 CALL SET_OUTFLOW_FLUXES(IJK, FIJK)
178 ENDIF
179
180
181
182
183 IF (FLUID_AT(JM_OF(IJK))) THEN
184 FIJK = JM_OF(IJK)
185 RVEL_G = V_G(FIJK)
186 IF(MMAX>0) RVEL_S(:MMAX) = V_S(FIJK,:MMAX)
187
188 CALL SET_OUTFLOW_MISC(BCV, IJK, FIJK)
189 CALL SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
190 IF (BC_TYPE(BCV)=='P_INFLOW') &
191 CALL SET_PINOUTFLOW(BCV, IJK, FIJK, RVEL_G, RVEL_S)
192
193 IF (ROP_G(IJK) > ZERO) THEN
194 V_G(IJK) = ROP_G(FIJK)*V_G(FIJK)/ROP_G(IJK)
195 ELSE
196 V_G(IJK) = ZERO
197 ENDIF
198 U_G(IJK) = U_G(FIJK)
199 W_G(IJK) = W_G(FIJK)
200
201 IF (MMAX > 0) THEN
202 WHERE (ROP_S(IJK,:MMAX) > ZERO)
203 V_S(IJK,:MMAX) = ROP_S(FIJK,:MMAX)*&
204 V_S(FIJK,:MMAX)/ROP_S(IJK,:MMAX)
205 ELSEWHERE
206 V_S(IJK,:MMAX) = ZERO
207 END WHERE
208 U_S(IJK,:MMAX) = U_S(FIJK,:MMAX)
209 W_S(IJK,:MMAX) = W_S(FIJK,:MMAX)
210 ENDIF
211
212 CALL SET_OUTFLOW_FLUXES(IJK, FIJK)
213 ENDIF
214
215
216
217
218 IF (FLUID_AT(JP_OF(IJK))) THEN
219 FIJK = JP_OF(IJK)
220 RVEL_G = -V_G(IJK)
221 IF (MMAX>0) RVEL_S(:MMAX) = -V_S(IJK,:MMAX)
222
223 CALL SET_OUTFLOW_MISC(BCV, IJK, FIJK)
224 CALL SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
225 IF (BC_TYPE(BCV)=='P_INFLOW') &
226 CALL SET_PINOUTFLOW(BCV, IJK, FIJK, RVEL_G, RVEL_S)
227
228 IF (V_G(IJK) == UNDEFINED) THEN
229 IF (ROP_G(IJK) > ZERO) THEN
230 V_G(IJK) = ROP_G(FIJK)*V_G(FIJK)/ROP_G(IJK)
231 ELSE
232 V_G(IJK) = ZERO
233 ENDIF
234 ENDIF
235 U_G(IJK) = U_G(FIJK)
236 W_G(IJK) = W_G(FIJK)
237
238 DO M = 1, MMAX
239 IF (V_S(IJK,M) == UNDEFINED) THEN
240 IF (ROP_S(IJK,M) > ZERO) THEN
241 V_S(IJK,M) = ROP_S(FIJK,M)*&
242 V_S(FIJK,M)/ROP_S(IJK,M)
243 ELSE
244 V_S(IJK,M) = ZERO
245 ENDIF
246 ENDIF
247 U_S(IJK,M) = U_S(FIJK,M)
248 W_S(IJK,M) = W_S(FIJK,M)
249 ENDDO
250
251 CALL SET_OUTFLOW_FLUXES(IJK, FIJK)
252 ENDIF
253
254
255
256
257 IF (FLUID_AT(KM_OF(IJK))) THEN
258 FIJK = KM_OF(IJK)
259 RVEL_G = W_G(FIJK)
260 IF (MMAX>0) RVEL_S(:MMAX) = W_S(FIJK,:MMAX)
261
262 CALL SET_OUTFLOW_MISC(BCV, IJK, FIJK)
263 CALL SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
264 IF (BC_TYPE(BCV)=='P_INFLOW') &
265 CALL SET_PINOUTFLOW(BCV, IJK, FIJK, RVEL_G, RVEL_S)
266
267 IF (ROP_G(IJK) > ZERO) THEN
268 W_G(IJK) = ROP_G(FIJK)*W_G(FIJK)/ROP_G(IJK)
269 ELSE
270 W_G(IJK) = ZERO
271 ENDIF
272 U_G(IJK) = U_G(FIJK)
273 V_G(IJK) = V_G(FIJK)
274
275 IF (MMAX > 0) THEN
276 WHERE (ROP_S(IJK,:MMAX) > ZERO)
277 W_S(IJK,:MMAX) = ROP_S(FIJK,:MMAX)*&
278 W_S(FIJK,:MMAX)/ROP_S(IJK,:MMAX)
279 ELSEWHERE
280 W_S(IJK,:MMAX) = ZERO
281 END WHERE
282 U_S(IJK,:MMAX) = U_S(FIJK,:MMAX)
283 V_S(IJK,:MMAX) = V_S(FIJK,:MMAX)
284 ENDIF
285
286 CALL SET_OUTFLOW_FLUXES(IJK, FIJK)
287 ENDIF
288
289
290
291
292 IF (FLUID_AT(KP_OF(IJK))) THEN
293 FIJK = KP_OF(IJK)
294 RVEL_G = -W_G(IJK)
295 IF (MMAX>0) RVEL_S(:MMAX) = -W_S(IJK,:MMAX)
296
297 CALL SET_OUTFLOW_MISC(BCV, IJK, FIJK)
298 CALL SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
299 IF (BC_TYPE(BCV)=='P_INFLOW') &
300 CALL SET_PINOUTFLOW(BCV, IJK, FIJK, RVEL_G, RVEL_S)
301
302 IF (W_G(IJK) == UNDEFINED) THEN
303 IF (ROP_G(IJK) > ZERO) THEN
304 W_G(IJK) = ROP_G(FIJK)*W_G(FIJK)/ROP_G(IJK)
305 ELSE
306 W_G(IJK) = ZERO
307 ENDIF
308 ENDIF
309 U_G(IJK) = U_G(FIJK)
310 V_G(IJK) = V_G(FIJK)
311
312 DO M = 1, MMAX
313 IF (W_S(IJK,M) == UNDEFINED) THEN
314 IF (ROP_S(IJK,M) > ZERO) THEN
315 W_S(IJK,M) = ROP_S(FIJK,M)*&
316 W_S(FIJK,M)/ROP_S(IJK,M)
317 ELSE
318 W_S(IJK,M) = ZERO
319 ENDIF
320 ENDIF
321 U_S(IJK,M) = U_S(FIJK,M)
322 V_S(IJK,M) = V_S(FIJK,M)
323 ENDDO
324
325 CALL SET_OUTFLOW_FLUXES(IJK, FIJK)
326 ENDIF
327
328 ENDDO
329 ENDDO
330 ENDDO
331
332 RETURN
333 END SUBROUTINE SET_OUTFLOW
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356 SUBROUTINE SET_OUTFLOW_MISC(BCV, IJK, FIJK)
357
358
359
360 use bc, only: bc_type
361 use run, only: kt_type_enum, ghd_2007
362 use fldvar, only: p_g, ro_g, T_g
363 use fldvar, only: p_s, p_star
364 use physprop, only: smax, mmax
365 use physprop, only: ro_g0, mw_mix_g
366 use eos, only: EOSG
367
368
369
370 use param1, only: undefined
371 IMPLICIT NONE
372
373
374
375
376 INTEGER, INTENT(IN) :: BCV
377
378 INTEGER, INTENT(IN) :: IJK
379
380 INTEGER, INTENT(IN) :: FIJK
381
382
383 IF (BC_TYPE(BCV) /= 'P_OUTFLOW' .AND. &
384 BC_TYPE(BCV) /= 'P_INFLOW') P_G(IJK) = P_G(FIJK)
385
386 MW_MIX_G(IJK) = MW_MIX_G(FIJK)
387
388
389
390 IF (RO_G0 == UNDEFINED) RO_G(IJK) = &
391 EOSG(MW_MIX_G(IJK),P_G(IJK),T_G(IJK))
392
393 IF (SMAX >0) THEN
394 P_STAR(IJK) = P_STAR(FIJK)
395 P_S(IJK,:SMAX) = P_S(FIJK,:SMAX)
396 ENDIF
397
398 IF (KT_TYPE_ENUM == GHD_2007 .AND. MMAX>0) &
399 P_S(IJK,MMAX) = P_S(FIJK,MMAX)
400
401 RETURN
402 END SUBROUTINE SET_OUTFLOW_MISC
403
404
405
406
407
408
409
410
411
412
413
414 SUBROUTINE SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
415
416
417
418 use bc, only: bc_type, bc_ep_g, bc_rop_s
419 use run, only: kt_type_enum, ghd_2007
420 use physprop, only: smax, mmax
421 use fldvar, only: rop_g, ro_g, ep_g
422 use fldvar, only: rop_s, ro_s, ep_s
423 use discretelement, only: discrete_element, des_mmax
424 use discretelement, only: des_rop_s, des_ro_s
425
426
427 use constant, only: pi
428 use param, only: dimension_m
429 use param1, only: undefined, zero, one
430 IMPLICIT NONE
431
432
433
434
435 INTEGER, INTENT(IN) :: BCV
436
437 INTEGER, INTENT(IN) :: IJK
438
439 INTEGER, INTENT(IN) :: FIJK
440
441
442
443
444
445
446
447 DOUBLE PRECISION, INTENT(IN) :: RVEL_G
448 DOUBLE PRECISION, INTENT(IN), DIMENSION(DIMENSION_M) :: RVEL_S
449
450
451
452
453 INTEGER :: M, mm
454
455 DOUBLE PRECISION :: EPs
456
457 DOUBLE PRECISION :: SUM_EPs
458
459 DOUBLE PRECISION :: SUM_ROPS
460
461
462
463 = ZERO
464 SUM_EPS = ZERO
465
466 DO M = 1, SMAX
467
468 IF(BC_TYPE(BCV) == 'P_INFLOW') THEN
469 ROP_S(IJK,M) = ROP_S(FIJK,M)
470 ELSE
471
472
473 IF (RVEL_S(M) >= ZERO) THEN
474
475 (IJK,M) = ROP_S(FIJK,M)
476 ELSE
477
478 (IJK,M) = ZERO
479 ENDIF
480 ENDIF
481
482
483
484 IF(BC_ROP_S(BCV,M)/=UNDEFINED) ROP_S(IJK,M)=BC_ROP_S(BCV,M)
485
486
487 = SUM_ROPS + ROP_S(IJK,M)
488 SUM_EPS = SUM_EPS + EP_S(IJK,M)
489 ENDDO
490
491 IF (KT_TYPE_ENUM == GHD_2007) ROP_S(IJK,MMAX) = SUM_ROPS
492
493
494
495
496 IF (DISCRETE_ELEMENT .AND. ALLOCATED(DES_ROP_S)) THEN
497 DO M = 1, DES_MMAX
498
499
500
501
502 (IJK,M) = DES_ROP_S(FIJK,M)
503 SUM_ROPS = SUM_ROPS + DES_ROP_S(IJK,M)
504 EPS = DES_ROP_S(IJK,M)/DES_RO_S(M)
505 SUM_EPS = SUM_EPS + EPS
506 ENDDO
507 ENDIF
508
509
510
511
512 IF (BC_EP_G(BCV) == UNDEFINED) EP_G(IJK) = ONE - SUM_EPS
513
514
515
516 (IJK) = RO_G(IJK)*EP_G(IJK)
517
518 RETURN
519 END SUBROUTINE SET_OUTFLOW_EP
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537 SUBROUTINE SET_OUTFLOW_FLUXES(IJK,FIJK)
538
539
540
541 use physprop, only: mmax
542 use run, only: kt_type_enum, ghd_2007
543 use run, only: added_mass
544 use mflux, only: flux_ge, flux_gn, flux_gt
545 use mflux, only: flux_se, flux_sn, flux_st
546 use mflux, only: flux_ne, flux_nn, flux_nt
547 use mflux, only: flux_gse, flux_gsn, flux_gst
548 use mflux, only: flux_sse, flux_ssn, flux_sst
549
550
551
552
553 INTEGER, INTENT(IN) :: IJK
554
555 INTEGER, INTENT(IN) :: FIJK
556
557
558 (IJK) = Flux_gE(FIJK)
559 Flux_gN(IJK) = Flux_gN(FIJK)
560 Flux_gT(IJK) = Flux_gT(FIJK)
561 IF (MMAX >0) THEN
562 Flux_sE(IJK,:MMAX) = Flux_sE(FIJK,:MMAX)
563 Flux_sN(IJK,:MMAX) = Flux_sN(FIJK,:MMAX)
564 Flux_sT(IJK,:MMAX) = Flux_sT(FIJK,:MMAX)
565 ENDIF
566
567 IF(ADDED_MASS) THEN
568 Flux_gSE(IJK) = Flux_gSE(FIJK)
569 Flux_gSN(IJK) = Flux_gSN(FIJK)
570 Flux_gST(IJK) = Flux_gST(FIJK)
571 IF (MMAX >0) THEN
572 Flux_sSE(IJK) = Flux_sSE(FIJK)
573 Flux_sSN(IJK) = Flux_sSN(FIJK)
574 Flux_sST(IJK) = Flux_sST(FIJK)
575 ENDIF
576 ENDIF
577
578 IF (KT_TYPE_ENUM == GHD_2007) THEN
579 Flux_nE(IJK) = Flux_nE(FIJK)
580 Flux_nN(IJK) = Flux_nN(FIJK)
581 Flux_nT(IJK) = Flux_nT(FIJK)
582 ENDIF
583
584 RETURN
585 END SUBROUTINE SET_OUTFLOW_FLUXES
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602 SUBROUTINE SET_PINOUTFLOW(BCV,IJK,FIJK,RVEL_G,RVEL_S)
603
604
605
606 use bc, only: bc_type
607 use bc, only: bc_t_g, bc_x_g
608 use bc, only: bc_scalar, bc_k_turb_g, bc_e_turb_g
609 use bc, only: bc_t_s, bc_x_s, bc_theta_m
610 use fldvar, only: t_g, t_s, theta_m
611 use fldvar, only: x_g, x_s, scalar
612 use fldvar, only: k_turb_g, e_turb_g
613
614 use fldvar, only: rop_s, ro_s, d_p
615 use run, only: kt_type_enum, ghd_2007
616 use run, only: k_epsilon
617 use physprop, only: nmax, smax, mmax
618 use scalars, only: nscalar, phase4scalar
619
620
621
622 use param, only: dimension_m
623 use param1, only: zero
624 use constant, only: pi
625 IMPLICIT NONE
626
627
628
629
630 INTEGER, INTENT(IN) :: BCV
631
632 INTEGER, INTENT(IN) :: IJK
633
634 INTEGER, INTENT(IN) :: FIJK
635
636
637 DOUBLE PRECISION, INTENT(IN) :: RVEL_G, RVEL_S(DIMENSION_M)
638
639
640
641
642 DOUBLE PRECISION :: Nm, NTot
643
644 INTEGER :: M, N, Ms
645
646
647
648 IF (RVEL_G < 0) THEN
649 (IJK) = BC_T_G(BCV)
650 IF (NMAX(0) > 0) &
651 X_G(IJK,:NMAX(0)) = BC_X_G(BCV,:NMAX(0))
652 IF (K_Epsilon) THEN
653 K_Turb_G(IJK) = BC_K_Turb_G(BCV)
654 E_Turb_G(IJK) = BC_E_Turb_G(BCV)
655 ENDIF
656 IF (NScalar > 0) THEN
657 DO N = 1,NScalar
658 Ms = Phase4Scalar(N)
659 IF (Ms == 0) &
660 Scalar(IJK, N) = BC_Scalar(BCV,N)
661 ENDDO
662 ENDIF
663 ELSE
664 (IJK) = T_G(FIJK)
665 IF (NMAX(0) > 0) &
666 X_G(IJK,:NMAX(0)) = X_G(FIJK,:NMAX(0))
667 IF(K_Epsilon) THEN
668 K_Turb_G(IJK) = K_Turb_G(FIJK)
669 E_Turb_G(IJK) = E_Turb_G(FIJK)
670 ENDIF
671 IF (NScalar >0) THEN
672 DO N = 1, NScalar
673 Ms = Phase4Scalar(N)
674 IF (Ms == 0) &
675 Scalar(IJK, N) = Scalar(FIJK, N)
676 ENDDO
677 ENDIF
678 ENDIF
679
680
681 DO M = 1, SMAX
682 IF (RVEL_S(M) < 0) THEN
683 (IJK,M) = BC_T_S(BCV,M)
684 THETA_M(IJK,M) = BC_THETA_M(BCV,M)
685 IF (NMAX(M) > 0) &
686 X_S(IJK,M,:NMAX(M)) = BC_X_S(BCV,M,:NMAX(M))
687 IF (NScalar > 0) THEN
688 DO N = 1,NScalar
689 Ms = Phase4Scalar(N)
690 IF (Ms == M) &
691 Scalar(IJK, N) = BC_Scalar(BCV,N)
692 ENDDO
693 ENDIF
694 ELSE
695 (IJK,M) = T_S(FIJK,M)
696 THETA_M(IJK,M) = THETA_M(FIJK,M)
697 IF (NMAX(M) > 0) &
698 X_S(IJK,M,:NMAX(M)) = X_S(FIJK,M,:NMAX(M))
699 IF (NScalar > 0) THEN
700 DO N = 1,NScalar
701 Ms = Phase4Scalar(N)
702 IF (Ms == M) &
703 Scalar(IJK, N) = Scalar(FIJK,N)
704 ENDDO
705 ENDIF
706 ENDIF
707 ENDDO
708
709
710
711 IF(KT_TYPE_ENUM == GHD_2007) THEN
712 nTOT = zero
713 DO M = 1, SMAX
714 nM = ROP_S(IJK,M)*6d0/ &
715 (PI*D_p(IJK,M)**3*RO_S(IJK,M))
716 nTOT = nTOT + nM
717 THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) + &
718 nM*THETA_M(IJK,M)
719 ENDDO
720 IF (NTOT > zero) &
721 THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) / nTOT
722 ENDIF
723
724 RETURN
725 END SUBROUTINE SET_PINOUTFLOW
726