File: /nfs/home/0/users/jenkins/mfix.git/model/cartesian_grid/CG_source_w_g.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21 SUBROUTINE CG_SOURCE_W_G(A_M, B_M)
22
23
24
25
26
27
28
29
30 USE param
31 USE param1
32 USE parallel
33 USE matrix
34 USE scales
35 USE constant
36 USE physprop
37 USE fldvar
38 USE visc_g
39 USE rxns
40 USE run
41 USE toleranc
42 USE geometry
43 USE indices
44 USE is
45 USE tau_g
46 USE bc
47 USE compar
48 USE sendrecv
49 USE ghdtheory
50 USE drag
51 USE fun_avg
52 USE functions
53
54
55
56 USE cutcell
57 USE quadric
58
59
60
61
62 IMPLICIT NONE
63
64
65
66
67
68
69
70
71
72 INTEGER IER
73
74
75 INTEGER I, J, K, IJK, IJKT, IMJK, IJKP, IMJKP,&
76 IJKE, IJKW, IJKTE, IM, IPJK
77
78
79 INTEGER M
80
81
82 DOUBLE PRECISION EPGA
83
84
85 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
86
87
88 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
89
90
91
92
93 INTEGER :: JM,IP,JP,IJMK,IJPK,IJKC,IJKN,IJKNE,IJKS,IJKSE,IPJMK,IJKM,KM,KP,IJMKP
94 INTEGER :: IJKTN,IJKWT,IJKST
95 DOUBLE PRECISION :: We,Ww,Wn,Ws,Wt,Wb
96 DOUBLE PRECISION :: B_NOC
97 DOUBLE PRECISION :: MU_GT_E,MU_GT_W,MU_GT_N,MU_GT_S,MU_GT_T,MU_GT_B,MU_GT_CUT
98 DOUBLE PRECISION :: WW_g
99 INTEGER :: BCV
100 CHARACTER(LEN=9) :: BCT
101
102 DOUBLE PRECISION F_vir, ROP_MA, U_se, Usw, Wse, Wsw, Wsn, Wss, Wst, Wsb, Usc,Vsc,Vsn,Vss
103
104 DOUBLE PRECISION :: W_F_Slip
105
106
107
108
109
110 IF(CG_SAFE_MODE(5)==1) RETURN
111
112
113 = 0
114 IF (.NOT.MOMENTUM_Z_EQ(0)) RETURN
115
116
117
118
119
120
121 DO IJK = ijkstart3, ijkend3
122 I = I_OF(IJK)
123 J = J_OF(IJK)
124 K = K_OF(IJK)
125 IJKT = TOP_OF(IJK)
126 EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
127 IF (IP_AT_T(IJK)) THEN
128
129
130
131
132 ELSE IF (EPGA <= DIL_EP_S) THEN
133
134
135
136
137 ELSEIF(INTERIOR_CELL_AT(IJK)) THEN
138
139 = BC_W_ID(IJK)
140
141 IF(BCV > 0 ) THEN
142 BCT = BC_TYPE(BCV)
143 ELSE
144 BCT = 'NONE'
145 ENDIF
146
147 SELECT CASE (BCT)
148 CASE ('CG_NSW')
149 NOC_WG = .TRUE.
150 WW_g = ZERO
151 MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJKT)*MU_GT(IJKT))/(VOL(IJK) + VOL(IJKT))
152
153 IF(.NOT.K_EPSILON) THEN
154 A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_W_CUT(IJK)/DELH_W(IJK)
155 ELSE
156 CALL Wall_Function(IJK,IJK,ONE/DELH_W(IJK),W_F_Slip)
157 A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_W_CUT(IJK)*(ONE-W_F_Slip)/DELH_W(IJK)
158 ENDIF
159
160 CASE ('CG_FSW')
161 NOC_WG = .FALSE.
162 WW_g = ZERO
163 CASE('CG_PSW')
164 IF(BC_HW_G(BCV)==UNDEFINED) THEN
165 = .TRUE.
166 WW_g = BC_WW_G(BCV)
167 MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJKT)*MU_GT(IJKT))/(VOL(IJK) + VOL(IJKT))
168 A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_W_CUT(IJK)/DELH_W(IJK)
169 B_M(IJK,M) = B_M(IJK,M) - MU_GT_CUT * WW_g * Area_W_CUT(IJK)/DELH_W(IJK)
170 ELSEIF(BC_HW_G(BCV)==ZERO) THEN
171 = .FALSE.
172 WW_g = ZERO
173 ELSE
174 = .FALSE.
175 WW_g = BC_WW_G(BCV)
176 MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJKT)*MU_GT(IJKT))/(VOL(IJK) + VOL(IJKT))
177 A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_W_CUT(IJK)*(BC_HW_G(BCV))
178 B_M(IJK,M) = B_M(IJK,M) - MU_GT_CUT * WW_g * Area_W_CUT(IJK)*(BC_HW_G(BCV))
179
180 ENDIF
181 CASE ('NONE')
182 NOC_WG = .FALSE.
183 END SELECT
184
185 IF(NOC_WG) THEN
186
187 IMJK = IM_OF(IJK)
188 IJMK = JM_OF(IJK)
189 IJKM = KM_OF(IJK)
190 IPJK = IP_OF(IJK)
191 IJPK = JP_OF(IJK)
192 IJKP = KP_OF(IJK)
193
194 We = Theta_We_bar(IJK) * W_g(IJK) + Theta_We(IJK) * W_g(IPJK)
195 Ww = Theta_We_bar(IMJK) * W_g(IMJK) + Theta_We(IMJK) * W_g(IJK)
196
197 Wn = Theta_Wn_bar(IJK) * W_g(IJK) + Theta_Wn(IJK) * W_g(IJPK)
198 Ws = Theta_Wn_bar(IJMK) * W_g(IJMK) + Theta_Wn(IJMK) * W_g(IJK)
199
200 Wt = Theta_Wt_bar(IJK) * W_g(IJK) + Theta_Wt(IJK) * W_g(IJKP)
201 Wb = Theta_Wt_bar(IJKM) * W_g(IJKM) + Theta_Wt(IJKM) * W_g(IJK)
202
203 IJKE = EAST_OF(IJK)
204
205 ijkt = top_of(ijk)
206
207 IF (WALL_AT(IJK)) THEN
208 IJKC = IJKT
209 ELSE
210 IJKC = IJK
211 ENDIF
212
213 IP = IP1(I)
214 IM = IM1(I)
215 IJKN = NORTH_OF(IJK)
216 IJKNE = EAST_OF(IJKN)
217
218 JM = JM1(J)
219 IPJMK = IP_OF(IJMK)
220 IJKS = SOUTH_OF(IJK)
221 IJKSE = EAST_OF(IJKS)
222
223 KP = KP1(K)
224 IJKT = TOP_OF(IJK)
225 IJKE = EAST_OF(IJK)
226 IJKP = KP_OF(IJK)
227 IJKTN = NORTH_OF(IJKT)
228 IJKTE = EAST_OF(IJKT)
229 IJKW = WEST_OF(IJK)
230 IJKWT = TOP_OF(IJKW)
231 IJKS = SOUTH_OF(IJK)
232 IJKST = TOP_OF(IJKS)
233
234 MU_GT_E = AVG_Z_H(AVG_X_H(MU_GT(IJKC),MU_GT(IJKE),I),&
235 AVG_X_H(MU_GT(IJKT),MU_GT(IJKTE),I),K)
236
237 MU_GT_W = AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJKC),IM),&
238 AVG_X_H(MU_GT(IJKWT),MU_GT(IJKT),IM),K)
239
240 MU_GT_N = AVG_Z_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),&
241 AVG_Y_H(MU_GT(IJKT),MU_GT(IJKTN),J),K)
242
243 MU_GT_S = AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),&
244 AVG_Y_H(MU_GT(IJKST),MU_GT(IJKT),JM),K)
245
246 MU_GT_T = MU_GT(IJKT)
247 MU_GT_B = MU_GT(IJKC)
248
249 B_NOC = MU_GT_E * Ayz_W(IJK) * (We-WW_g) * NOC_W_E(IJK) &
250 - MU_GT_W * Ayz_W(IMJK) * (Ww-WW_g) * NOC_W_E(IMJK) &
251 + MU_GT_N * Axz_W(IJK) * (Wn-WW_g) * NOC_W_N(IJK) &
252 - MU_GT_S * Axz_W(IJMK) * (Ws-WW_g) * NOC_W_N(IJMK) &
253 + MU_GT_T * Axy_W(IJK) * (Wt-WW_g) * NOC_W_T(IJK) *2.0d0&
254 - MU_GT_B * Axy_W(IJKM) * (Wb-WW_g) * NOC_W_T(IJKM) *2.0D0
255
256 B_M(IJK,M) = B_M(IJK,M) + B_NOC
257 ENDIF
258
259 IF(CUT_W_TREATMENT_AT(IJK)) THEN
260
261
262
263 = ZERO
264 IF(Added_Mass) THEN
265
266 F_vir = ( (W_s(IJK,M_AM) - W_sO(IJK,M_AM)) )*ODT*VOL_W(IJK)
267
268 I = I_OF(IJK)
269 J = J_OF(IJK)
270 K = K_OF(IJK)
271
272 IM = I - 1
273 JM = J - 1
274 KM = K - 1
275
276 IP = I + 1
277 JP = J + 1
278 KP = K + 1
279
280 IMJK = FUNIJK(IM,J,K)
281 IJMK = FUNIJK(I,JM,K)
282 IPJK = FUNIJK(IP,J,K)
283 IJPK = FUNIJK(I,JP,K)
284 IJKP = FUNIJK(I,J,KP)
285 IJKM = FUNIJK(I,J,KM)
286
287 IMJKP = KP_OF(IMJK)
288 IJMKP = KP_OF(IJMK)
289
290 IJKE = EAST_OF(IJK)
291
292
293
294
295 = Theta_We_bar(IJK) * W_s(IJK,M_AM) + Theta_We(IJK) * W_s(IPJK,M_AM)
296 Wsw = Theta_We_bar(IMJK) * W_s(IMJK,M_AM) + Theta_We(IMJK) * W_s(IJK,M_AM)
297
298 U_se = Theta_W_te(IJK) * U_s(IJK,M_AM) + Theta_W_be(IJK) * U_s(IJKP,M_AM)
299 Usw = Theta_W_te(IMJK) * U_s(IMJK,M_AM) + Theta_W_be(IMJK) * U_s(IMJKP,M_AM)
300
301 Usc = (DELX_we(IJK) * Usw + DELX_ww(IJK) * U_se) / (DELX_we(IJK) + DELX_ww(IJK))
302
303
304 Wsn = Theta_Wn_bar(IJK) * W_s(IJK,M_AM) + Theta_Wn(IJK) * W_s(IJPK,M_AM)
305 Wss = Theta_Wn_bar(IJMK) * W_s(IJMK,M_AM) + Theta_Wn(IJMK) * W_s(IJK,M_AM)
306
307 Vsn = Theta_W_tn(IJK) * V_s(IJK,M_AM) + Theta_W_bn(IJK) * V_s(IJKP,M_AM)
308 Vss = Theta_W_tn(IJMK) * V_s(IJMK,M_AM) + Theta_W_bn(IJMK) * V_s(IJMKP,M_AM)
309
310 Vsc = (DELY_wn(IJK) * Vss + DELY_ws(IJK) * Vsn) / (DELY_wn(IJK) + DELY_ws(IJK))
311
312
313 Wst = Theta_Wt_bar(IJK) * W_s(IJK,M_AM) + Theta_Wt(IJK) * W_s(IJKP,M_AM)
314 Wsb = Theta_Wt_bar(IMJK) * W_s(IMJK,M_AM) + Theta_Wt(IMJK) * W_s(IMJKP,M_AM)
315
316
317
318
319 = F_vir + Usc * (Wse - Wsw)*AYZ(IJK) + &
320 Vsc * (Wsn - Wss)*AXZ(IJK) + &
321 W_s(IJK,M_AM)*(Wst - Wsb)*AXY(IJK)
322
323
324 ROP_MA = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M_AM) + VOL(IJKT)*ROP_g(IJKT)*EP_s(IJKT,M_AM))/(VOL(IJK) + VOL(IJKT))
325
326 F_vir = F_vir * Cv * ROP_MA
327
328 B_M(IJK,M) = B_M(IJK,M) - F_vir
329
330 ENDIF
331
332
333
334 ENDIF
335
336 ENDIF
337 END DO
338
339 RETURN
340 END SUBROUTINE CG_SOURCE_W_G
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363 SUBROUTINE CG_SOURCE_W_G_BC(A_M, B_M)
364
365
366
367
368
369
370
371
372 USE param
373 USE param1
374 USE parallel
375 USE matrix
376 USE scales
377 USE constant
378 USE physprop
379 USE fldvar
380 USE visc_g
381 USE rxns
382 USE run
383 USE toleranc
384 USE geometry
385 USE indices
386 USE is
387 USE tau_g
388 USE bc
389 USE output
390 USE compar
391 USE fun_avg
392 USE functions
393
394
395
396
397 USE cutcell
398 USE quadric
399
400
401
402
403
404 IMPLICIT NONE
405
406
407
408
409
410
411
412
413 INTEGER IJK, IJKB
414
415
416 INTEGER M
417
418
419 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
420
421
422 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
423
424
425
426
427
428 INTEGER :: BCV
429 CHARACTER(LEN=9) :: BCT
430
431
432
433 IF(CG_SAFE_MODE(5)==1) RETURN
434
435
436 = 0
437
438 DO IJK = ijkstart3, ijkend3
439
440 BCV = BC_W_ID(IJK)
441
442 IF(BCV > 0 ) THEN
443 BCT = BC_TYPE(BCV)
444 ELSE
445 BCT = 'NONE'
446 ENDIF
447
448 SELECT CASE (BCT)
449
450 CASE ('CG_NSW')
451
452 IF(WALL_W_AT(IJK)) THEN
453
454 A_M(IJK,E,M) = ZERO
455 A_M(IJK,W,M) = ZERO
456 A_M(IJK,N,M) = ZERO
457 A_M(IJK,S,M) = ZERO
458 A_M(IJK,T,M) = ZERO
459 A_M(IJK,B,M) = ZERO
460 A_M(IJK,0,M) = -ONE
461
462 B_M(IJK,M) = ZERO
463
464 ENDIF
465
466 CASE ('CG_FSW')
467
468 IF(WALL_W_AT(IJK)) THEN
469
470 A_M(IJK,E,M) = ZERO
471 A_M(IJK,W,M) = ZERO
472 A_M(IJK,N,M) = ZERO
473 A_M(IJK,S,M) = ZERO
474 A_M(IJK,T,M) = ZERO
475 A_M(IJK,B,M) = ZERO
476 A_M(IJK,0,M) = -ONE
477
478
479
480 (IJK,M) = ZERO
481
482 IF(DABS(NORMAL_W(IJK,3))/=ONE) THEN
483
484 IF (W_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
485 A_M(IJK,E,M) = ONE
486 ELSEIF (W_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
487 A_M(IJK,W,M) = ONE
488 ELSEIF (W_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
489 A_M(IJK,N,M) = ONE
490 ELSEIF (W_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
491 A_M(IJK,S,M) = ONE
492 ELSEIF (W_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
493 A_M(IJK,T,M) = ONE
494 ELSEIF (W_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
495 A_M(IJK,B,M) = ONE
496 ENDIF
497
498 ENDIF
499
500 ENDIF
501
502 CASE ('CG_PSW')
503
504 IF(WALL_W_AT(IJK)) THEN
505
506 A_M(IJK,E,M) = ZERO
507 A_M(IJK,W,M) = ZERO
508 A_M(IJK,N,M) = ZERO
509 A_M(IJK,S,M) = ZERO
510 A_M(IJK,T,M) = ZERO
511 A_M(IJK,B,M) = ZERO
512 A_M(IJK,0,M) = -ONE
513
514
515 IF(BC_HW_G(BCV)==UNDEFINED) THEN
516 (IJK,M) = -BC_WW_G(BCV)
517 ELSEIF(BC_HW_G(BCV)==ZERO) THEN
518 (IJK,M) = ZERO
519
520 IF(DABS(NORMAL_W(IJK,3))/=ONE) THEN
521
522 IF (W_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
523 A_M(IJK,E,M) = ONE
524 ELSEIF (W_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
525 A_M(IJK,W,M) = ONE
526 ELSEIF (W_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
527 A_M(IJK,N,M) = ONE
528 ELSEIF (W_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
529 A_M(IJK,S,M) = ONE
530 ELSEIF (W_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
531 A_M(IJK,T,M) = ONE
532 ELSEIF (W_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
533 A_M(IJK,B,M) = ONE
534 ENDIF
535
536 ENDIF
537
538 ELSE
539
540
541
542
543
544 ENDIF
545
546 ENDIF
547
548
549 CASE ('CG_MI')
550
551 A_M(IJK,E,M) = ZERO
552 A_M(IJK,W,M) = ZERO
553 A_M(IJK,N,M) = ZERO
554 A_M(IJK,S,M) = ZERO
555 A_M(IJK,T,M) = ZERO
556 A_M(IJK,B,M) = ZERO
557 A_M(IJK,0,M) = -ONE
558
559 IF(BC_W_g(BCV)/=UNDEFINED) THEN
560 B_M(IJK,M) = - BC_W_g(BCV)
561 ELSE
562 B_M(IJK,M) = - BC_VELMAG_g(BCV)*NORMAL_W(IJK,3)
563 ENDIF
564
565
566 IJKB = BOTTOM_OF(IJK)
567 IF(FLUID_AT(IJKB)) THEN
568
569 A_M(IJKB,E,M) = ZERO
570 A_M(IJKB,W,M) = ZERO
571 A_M(IJKB,N,M) = ZERO
572 A_M(IJKB,S,M) = ZERO
573 A_M(IJKB,T,M) = ZERO
574 A_M(IJKB,B,M) = ZERO
575 A_M(IJKB,0,M) = -ONE
576
577 IF(BC_W_g(BCV)/=UNDEFINED) THEN
578 B_M(IJKB,M) = - BC_W_g(BCV)
579 ELSE
580 B_M(IJKB,M) = - BC_VELMAG_g(BCV)*NORMAL_W(IJK,3)
581 ENDIF
582
583
584 ENDIF
585
586 CASE ('CG_PO')
587
588 A_M(IJK,E,M) = ZERO
589 A_M(IJK,W,M) = ZERO
590 A_M(IJK,N,M) = ZERO
591 A_M(IJK,S,M) = ZERO
592 A_M(IJK,T,M) = ZERO
593 A_M(IJK,B,M) = ZERO
594 A_M(IJK,0,M) = -ONE
595 B_M(IJK,M) = ZERO
596
597 IJKB = BOTTOM_OF(IJK)
598 IF(FLUID_AT(IJKB)) THEN
599
600 A_M(IJK,B,M) = ONE
601 A_M(IJK,0,M) = -ONE
602
603 ENDIF
604
605 END SELECT
606
607 BCV = BC_ID(IJK)
608
609 IF(BCV > 0 ) THEN
610 BCT = BC_TYPE(BCV)
611 ELSE
612 BCT = 'NONE'
613 ENDIF
614
615 SELECT CASE (BCT)
616
617 CASE ('CG_MI')
618
619 A_M(IJK,E,M) = ZERO
620 A_M(IJK,W,M) = ZERO
621 A_M(IJK,N,M) = ZERO
622 A_M(IJK,S,M) = ZERO
623 A_M(IJK,T,M) = ZERO
624 A_M(IJK,B,M) = ZERO
625 A_M(IJK,0,M) = -ONE
626
627 IF(BC_W_g(BCV)/=UNDEFINED) THEN
628 B_M(IJK,M) = - BC_W_g(BCV)
629 ELSE
630 B_M(IJK,M) = - BC_VELMAG_g(BCV)*NORMAL_S(IJK,3)
631 ENDIF
632
633
634 IJKB = BOTTOM_OF(IJK)
635 IF(FLUID_AT(IJKB)) THEN
636
637 A_M(IJKB,E,M) = ZERO
638 A_M(IJKB,W,M) = ZERO
639 A_M(IJKB,N,M) = ZERO
640 A_M(IJKB,S,M) = ZERO
641 A_M(IJKB,T,M) = ZERO
642 A_M(IJKB,B,M) = ZERO
643 A_M(IJKB,0,M) = -ONE
644
645 IF(BC_W_g(BCV)/=UNDEFINED) THEN
646 B_M(IJKB,M) = - BC_W_g(BCV)
647 ELSE
648 B_M(IJKB,M) = - BC_VELMAG_g(BCV)*NORMAL_S(IJK,3)
649 ENDIF
650
651
652 ENDIF
653
654 CASE ('CG_PO')
655
656 A_M(IJK,E,M) = ZERO
657 A_M(IJK,W,M) = ZERO
658 A_M(IJK,N,M) = ZERO
659 A_M(IJK,S,M) = ZERO
660 A_M(IJK,T,M) = ZERO
661 A_M(IJK,B,M) = ZERO
662 A_M(IJK,0,M) = -ONE
663 B_M(IJK,M) = ZERO
664
665 IJKB = BOTTOM_OF(IJK)
666 IF(FLUID_AT(IJKB)) THEN
667
668 A_M(IJK,B,M) = ONE
669 A_M(IJK,0,M) = -ONE
670
671 ENDIF
672
673 END SELECT
674
675
676 ENDDO
677
678 RETURN
679
680
681
682
683
684 END SUBROUTINE CG_SOURCE_W_G_BC
685
686
687
688
689
690