File: /nfs/home/0/users/jenkins/mfix.git/model/source_u_g.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 SUBROUTINE SOURCE_U_G(A_M, B_M, IER)
25
26
27
28
29 USE bc
30 USE bodyforce
31 USE compar
32 USE constant
33 USE cutcell
34 USE drag
35 USE fldvar
36 USE fun_avg
37 USE functions
38 USE geometry
39 USE ghdtheory
40 USE indices
41 USE is
42 USE matrix
43 USE mms
44 USE parallel
45 USE param
46 USE param1
47 USE physprop
48 USE quadric
49 USE run
50 USE rxns
51 USE scales
52 USE sendrecv
53 USE tau_g
54 USE toleranc
55 USE visc_g
56
57 IMPLICIT NONE
58
59
60
61
62 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
63
64 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
65
66 INTEGER, INTENT(INOUT) :: IER
67
68
69
70
71 INTEGER :: I, J, K, IJK, IJKE, IPJK, IJKM, &
72 IPJKM, IMJK, IJMK, IPJMK, IJPK, IJKP
73
74 INTEGER :: M, L, MM
75
76 INTEGER :: ISV
77
78 DOUBLE PRECISION :: PgE
79
80 DOUBLE PRECISION :: EPGA
81
82 DOUBLE PRECISION :: ROPGA, ROGA
83
84 DOUBLE PRECISION :: MUGA
85
86 DOUBLE PRECISION :: EPMUGA
87
88 DOUBLE PRECISION :: Wge
89
90 DOUBLE PRECISION :: Sdp
91
92 DOUBLE PRECISION :: V0, Vpm, Vmt, Vbf, Vcf, Vtza
93
94 DOUBLE PRECISION :: Ghd_drag, avgRop
95
96 DOUBLE PRECISION :: HYS_drag, avgDrag
97
98 DOUBLE PRECISION :: ROP_MA, U_se, Usw, Vsw, Vse, Usn,&
99 Uss, Wsb, Wst, Wse, Usb, Ust
100 DOUBLE PRECISION :: F_vir
101
102
103
104 = 0
105
106 IF (.NOT.MOMENTUM_X_EQ(0)) RETURN
107
108
109
110
111
112
113
114
115
116
117 DO IJK = ijkstart3, ijkend3
118 I = I_OF(IJK)
119 J = J_OF(IJK)
120 K = K_OF(IJK)
121 IJKE = EAST_OF(IJK)
122 IJKM = KM_OF(IJK)
123 IPJK = IP_OF(IJK)
124 IMJK = IM_OF(IJK)
125 IPJKM = IP_OF(IJKM)
126 IJMK = JM_OF(IJK)
127 IPJMK = IP_OF(IJMK)
128 IJPK = JP_OF(IJK)
129 IJKP = KP_OF(IJK)
130
131 EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
132
133
134 IF (IP_AT_E(IJK)) THEN
135 A_M(IJK,E,M) = ZERO
136 A_M(IJK,W,M) = ZERO
137 A_M(IJK,N,M) = ZERO
138 A_M(IJK,S,M) = ZERO
139 A_M(IJK,T,M) = ZERO
140 A_M(IJK,B,M) = ZERO
141 A_M(IJK,0,M) = -ONE
142 B_M(IJK,M) = ZERO
143
144
145 ELSEIF (EPGA <= DIL_EP_S) THEN
146 A_M(IJK,E,M) = ZERO
147 A_M(IJK,W,M) = ZERO
148 A_M(IJK,N,M) = ZERO
149 A_M(IJK,S,M) = ZERO
150 A_M(IJK,T,M) = ZERO
151 A_M(IJK,B,M) = ZERO
152 A_M(IJK,0,M) = -ONE
153 B_M(IJK,M) = ZERO
154
155
156 IF (EP_G(WEST_OF(IJK)) > DIL_EP_S) THEN
157 A_M(IJK,W,M) = ONE
158 ELSE IF (EP_G(EAST_OF(IJK)) > DIL_EP_S) THEN
159 A_M(IJK,E,M) = ONE
160 ELSE
161 B_M(IJK,M) = -U_G(IJK)
162 ENDIF
163
164
165 ELSEIF (BLOCKED_U_CELL_AT(IJK)) THEN
166 A_M(IJK,E,M) = ZERO
167 A_M(IJK,W,M) = ZERO
168 A_M(IJK,N,M) = ZERO
169 A_M(IJK,S,M) = ZERO
170 A_M(IJK,T,M) = ZERO
171 A_M(IJK,B,M) = ZERO
172 A_M(IJK,0,M) = -ONE
173 B_M(IJK,M) = ZERO
174
175
176 ELSE
177
178
179
180 = P_G(IJKE)
181 IF (CYCLIC_X_PD) THEN
182 IF (IMAP(I_OF(IJK)).EQ.IMAX1) PGE = P_G(IJKE) - DELP_X
183 ENDIF
184 IF (MODEL_B) THEN
185 IF(.NOT.CUT_U_TREATMENT_AT(IJK)) THEN
186 SDP = -P_SCALE*(PGE - P_G(IJK))*AYZ(IJK)
187 ELSE
188 SDP = -P_SCALE*(PGE * A_UPG_E(IJK) - P_G(IJK) * A_UPG_W(IJK) )
189 ENDIF
190 ELSE
191 IF(.NOT.CUT_U_TREATMENT_AT(IJK)) THEN
192 SDP = -P_SCALE*EPGA*(PGE - P_G(IJK))*AYZ(IJK)
193 ELSE
194 SDP = -P_SCALE*EPGA*(PGE * A_UPG_E(IJK) - P_G(IJK) * A_UPG_W(IJK) )
195 ENDIF
196 ENDIF
197
198 IF(.NOT.CUT_U_TREATMENT_AT(IJK)) THEN
199
200 = HALF * (VOL(IJK)*ROP_G(IJK) + &
201 VOL(IPJK)*ROP_G(IJKE))/VOL_U(IJK)
202 ROGA = HALF * (VOL(IJK)*RO_G(IJK) + &
203 VOL(IPJK)*RO_G(IJKE))/VOL_U(IJK)
204
205 = HALF * (VOL(IJK)*ROP_GO(IJK) + &
206 VOL(IPJK)*ROP_GO(IJKE))*ODT/VOL_U(IJK)
207
208 IF(Added_Mass) THEN
209 ROP_MA = AVG_X(ROP_g(IJK)*EP_s(IJK,M_AM),&
210 ROP_g(IJKE)*EP_s(IJKE,M_AM),I)
211 V0 = V0 + Cv * ROP_MA * ODT
212 ENDIF
213 ELSE
214
215 = (VOL(IJK)*ROP_G(IJK) + &
216 VOL(IPJK)*ROP_G(IJKE))/(VOL(IJK) + VOL(IPJK))
217 ROGA = (VOL(IJK)*RO_G(IJK) + &
218 VOL(IPJK)*RO_G(IJKE) )/(VOL(IJK) + VOL(IPJK))
219
220 = (VOL(IJK)*ROP_GO(IJK) + VOL(IPJK)*ROP_GO(IJKE))*&
221 ODT/(VOL(IJK) + VOL(IPJK))
222
223 IF(Added_Mass) THEN
224 ROP_MA = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M_AM) + &
225 VOL(IPJK)*ROP_g(IJKE)*EP_s(IJKE,M_AM) )/&
226 (VOL(IJK) + VOL(IPJK))
227 V0 = V0 + Cv * ROP_MA * ODT
228 ENDIF
229 ENDIF
230
231
232
233 = ZERO
234 IF(Added_Mass.AND.(.NOT.CUT_U_TREATMENT_AT(IJK))) THEN
235 F_vir = ( (U_s(IJK,M_AM) - U_sO(IJK,M_AM)) )*ODT*VOL_U(IJK)
236
237
238 = AVG_X_E(U_S(IMJK,M_AM),U_s(IJK,M_AM),I)
239 U_se = AVG_X_E(U_s(IJK,M_AM),U_s(IPJK,M_AM),IP1(I))
240 Vsw = AVG_Y_N(V_s(IJMK,M_AM),V_s(IJK,M_AM))
241 Vse = AVG_Y_N(V_s(IPJMK,M_AM),V_s(IPJK,M_AM))
242 Uss = AVG_Y(U_s(IJMK,M_AM),U_s(IJK,M_AM),JM1(J))
243 Usn = AVG_Y(U_s(IJK,M_AM),U_s(IJPK,M_AM),J)
244 IF(DO_K) THEN
245 Wsb = AVG_Z_T(W_s(IJKM,M_AM),W_s(IJK,M_AM))
246 Wst = AVG_Z_T(W_s(IPJKM,M_AM),W_s(IPJK,M_AM))
247 Wse = AVG_X(Wsb,Wst,I)
248 Usb = AVG_Z(U_s(IJKM,M_AM),U_s(IJK,M_AM),KM1(K))
249 Ust = AVG_Z(U_s(IJK,M_AM),U_s(IJKP,M_AM),K)
250 F_vir = F_vir + Wse*OX_E(I) * (Ust - Usb) *AXY(IJK)
251
252 IF (CYLINDRICAL) F_vir = F_vir - Wse**2*OX_E(I)
253 ENDIF
254
255 = F_vir + U_s(IJK,M_AM)*(U_se - Usw)*AYZ(IJK) + &
256 AVG_X(Vsw,Vse,I) * (Usn - Uss)*AXZ(IJK)
257 F_vir = F_vir * Cv * ROP_MA
258 ENDIF
259
260
261 IF (SIP_AT_E(IJK)) THEN
262 ISV = IS_ID_AT_E(IJK)
263 MUGA = AVG_X(MU_G(IJK),MU_G(IJKE),I)
264 VPM = MUGA/IS_PC(ISV,1)
265 IF (IS_PC(ISV,2) /= ZERO) VPM = VPM + &
266 HALF*IS_PC(ISV,2)*ROPGA*ABS(U_G(IJK))
267 ELSE
268 VPM = ZERO
269 ENDIF
270
271
272 IF(.NOT.CUT_U_TREATMENT_AT(IJK)) THEN
273 VMT = HALF * (VOL(IJK)*SUM_R_G(IJK) + &
274 VOL(IPJK)*SUM_R_G(IJKE))/VOL_U(IJK)
275 ELSE
276 VMT = (VOL(IJK)*SUM_R_G(IJK) + VOL(IPJK)*SUM_R_G(IJKE))/&
277 (VOL(IJK) + VOL(IPJK))
278 ENDIF
279
280
281 IF (MODEL_B) THEN
282 VBF = ROGA*BFX_G(IJK)
283 ELSE
284 = ROPGA*BFX_G(IJK)
285 ENDIF
286
287
288 = ZERO
289 IF (KT_TYPE_ENUM .EQ. GHD_2007) THEN
290 DO L = 1,SMAX
291 avgRop = AVG_X(ROP_S(IJK,L),ROP_S(IJKE,L),I)
292 if(avgRop > ZERO) Ghd_drag = Ghd_drag +&
293 AVG_X(F_GS(IJK,L),F_GS(IJKE,L),I) * JoiX(IJK,L) / avgRop
294 ENDDO
295 ENDIF
296
297
298 = ZERO
299 HYS_drag = ZERO
300 IF (DRAG_TYPE_ENUM .EQ. HYS .AND. KT_TYPE_ENUM .NE. GHD_2007) THEN
301 DO MM=1,MMAX
302 DO L = 1,MMAX
303 IF (L /= MM) THEN
304 avgDrag = AVG_X(beta_ij(IJK,MM,L),beta_ij(IJKE,MM,L),I)
305 HYS_drag = HYS_drag + avgDrag * (U_g(ijk) - U_s(IJK,L))
306 ENDIF
307 ENDDO
308 ENDDO
309 ENDIF
310
311
312 IF (CYLINDRICAL) THEN
313
314 = AVG_X(HALF*(W_G(IJK)+W_G(IJKM)),&
315 HALF*(W_G(IPJK)+W_G(IPJKM)),I)
316 VCF = ROPGA*WGE**2*OX_E(I)
317
318 IF(Added_Mass) VCF = VCF + Cv*ROP_MA*WGE**2*OX_E(I)
319
320
321 = AVG_X(MU_GT(IJK),MU_GT(IJKE),I)
322 VTZA = 2.d0*EPMUGA*OX_E(I)*OX_E(I)
323 ELSE
324 VCF = ZERO
325 VTZA = ZERO
326 ENDIF
327
328
329 (IJK,0,M) = -(A_M(IJK,E,M)+A_M(IJK,W,M)+&
330 A_M(IJK,N,M)+A_M(IJK,S,M)+A_M(IJK,T,M)+A_M(IJK,B,M)+&
331 (V0+VPM+ZMAX(VMT)+VTZA)*VOL_U(IJK))
332 B_M(IJK,M) = B_M(IJK,M) -(SDP + TAU_U_G(IJK) + &
333 ( (V0+ZMAX((-VMT)))*U_GO(IJK) + VBF + &
334 VCF + Ghd_drag + HYS_drag)*VOL_U(IJK) )
335
336 (IJK,M) = B_M(IJK,M) - F_vir
337
338 IF(USE_MMS) B_M(IJK,M) = &
339 B_M(IJK,M) - MMS_U_G_SRC(IJK)*VOL_U(IJK)
340 ENDIF
341 ENDDO
342
343
344
345 IF(CARTESIAN_GRID) CALL CG_SOURCE_U_G(A_M, B_M, IER)
346
347 CALL SOURCE_U_G_BC (A_M, B_M)
348
349 IF(CARTESIAN_GRID) CALL CG_SOURCE_U_G_BC(A_M, B_M, IER)
350
351 RETURN
352 END SUBROUTINE SOURCE_U_G
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373 SUBROUTINE SOURCE_U_G_BC(A_M, B_M)
374
375
376
377
378 USE param
379 USE param1
380 USE parallel
381 USE matrix
382 USE scales
383 USE constant
384 USE physprop
385 USE fldvar
386 USE visc_g
387 USE rxns
388 USE run
389 USE toleranc
390 USE geometry
391 USE indices
392 USE is
393 USE tau_g
394 USE bc
395 USE output
396 USE compar
397 USE fun_avg
398 USE functions
399 IMPLICIT NONE
400
401
402
403
404 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
405
406 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
407
408
409
410
411 INTEGER :: L
412
413 INTEGER :: I, J, K, IM, I1, I2, J1, J2, K1, K2, IJK,&
414 JM, KM, IJKW, IMJK, IP, IPJK
415
416 INTEGER :: M
417
418 DOUBLE PRECISION :: W_F_Slip
419
420
421
422 = 0
423
424
425
426
427
428
429
430
431
432
433
434 IF (DO_K) THEN
435
436 = 1
437 DO J1 = jmin3,jmax3
438 DO I1 = imin3, imax3
439 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
440 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE
441 = FUNIJK(I1,J1,K1)
442 IF (NS_WALL_AT(IJK)) THEN
443
444
445 (IJK,E,M) = ZERO
446 A_M(IJK,W,M) = ZERO
447 A_M(IJK,N,M) = ZERO
448 A_M(IJK,S,M) = ZERO
449 A_M(IJK,T,M) = -ONE
450 A_M(IJK,B,M) = ZERO
451 A_M(IJK,0,M) = -ONE
452 B_M(IJK,M) = ZERO
453 ELSEIF (FS_WALL_AT(IJK)) THEN
454
455
456 (IJK,E,M) = ZERO
457 A_M(IJK,W,M) = ZERO
458 A_M(IJK,N,M) = ZERO
459 A_M(IJK,S,M) = ZERO
460 A_M(IJK,T,M) = ONE
461 A_M(IJK,B,M) = ZERO
462 A_M(IJK,0,M) = -ONE
463 B_M(IJK,M) = ZERO
464 ENDIF
465 ENDDO
466 ENDDO
467
468
469 = KMAX2
470 DO J1 = jmin3,jmax3
471 DO I1 = imin3, imax3
472 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
473 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE
474 = FUNIJK(I1,J1,K1)
475 IF (NS_WALL_AT(IJK)) THEN
476 A_M(IJK,E,M) = ZERO
477 A_M(IJK,W,M) = ZERO
478 A_M(IJK,N,M) = ZERO
479 A_M(IJK,S,M) = ZERO
480 A_M(IJK,T,M) = ZERO
481 A_M(IJK,B,M) = -ONE
482 A_M(IJK,0,M) = -ONE
483 B_M(IJK,M) = ZERO
484 ELSEIF (FS_WALL_AT(IJK)) THEN
485 A_M(IJK,E,M) = ZERO
486 A_M(IJK,W,M) = ZERO
487 A_M(IJK,N,M) = ZERO
488 A_M(IJK,S,M) = ZERO
489 A_M(IJK,T,M) = ZERO
490 A_M(IJK,B,M) = ONE
491 A_M(IJK,0,M) = -ONE
492 B_M(IJK,M) = ZERO
493 ENDIF
494 ENDDO
495 ENDDO
496 ENDIF
497
498
499 = 1
500 DO K1 = kmin3, kmax3
501 DO I1 = imin3, imax3
502 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
503 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE
504 = FUNIJK(I1,J1,K1)
505 IF (NS_WALL_AT(IJK)) THEN
506 A_M(IJK,E,M) = ZERO
507 A_M(IJK,W,M) = ZERO
508 A_M(IJK,N,M) = -ONE
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 B_M(IJK,M) = ZERO
514 ELSEIF (FS_WALL_AT(IJK)) THEN
515 A_M(IJK,E,M) = ZERO
516 A_M(IJK,W,M) = ZERO
517 A_M(IJK,N,M) = ONE
518 A_M(IJK,S,M) = ZERO
519 A_M(IJK,T,M) = ZERO
520 A_M(IJK,B,M) = ZERO
521 A_M(IJK,0,M) = -ONE
522 B_M(IJK,M) = ZERO
523 ENDIF
524 ENDDO
525 ENDDO
526
527
528 = JMAX2
529 DO K1 = kmin3, kmax3
530 DO I1 = imin3, imax3
531 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
532 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE
533 = FUNIJK(I1,J1,K1)
534 IF (NS_WALL_AT(IJK)) THEN
535 A_M(IJK,E,M) = ZERO
536 A_M(IJK,W,M) = ZERO
537 A_M(IJK,N,M) = ZERO
538 A_M(IJK,S,M) = -ONE
539 A_M(IJK,T,M) = ZERO
540 A_M(IJK,B,M) = ZERO
541 A_M(IJK,0,M) = -ONE
542 B_M(IJK,M) = ZERO
543 ELSEIF (FS_WALL_AT(IJK)) THEN
544 A_M(IJK,E,M) = ZERO
545 A_M(IJK,W,M) = ZERO
546 A_M(IJK,N,M) = ZERO
547 A_M(IJK,S,M) = ONE
548 A_M(IJK,T,M) = ZERO
549 A_M(IJK,B,M) = ZERO
550 A_M(IJK,0,M) = -ONE
551 B_M(IJK,M) = ZERO
552 ENDIF
553 ENDDO
554 ENDDO
555
556
557
558
559
560 DO L = 1, DIMENSION_BC
561 IF (BC_DEFINED(L)) THEN
562
563
564
565 IF (BC_TYPE(L) == 'NO_SLIP_WALL' .AND. .NOT. K_Epsilon) THEN
566 I1 = BC_I_W(L)
567 I2 = BC_I_E(L)
568 J1 = BC_J_S(L)
569 J2 = BC_J_N(L)
570 K1 = BC_K_B(L)
571 K2 = BC_K_T(L)
572 DO K = K1, K2
573 DO J = J1, J2
574 DO I = I1, I2
575 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
576 IF (DEAD_CELL_AT(I,J,K)) CYCLE
577 = FUNIJK(I,J,K)
578 IF (.NOT.WALL_AT(IJK)) CYCLE
579 (IJK,E,M) = ZERO
580 A_M(IJK,W,M) = ZERO
581 A_M(IJK,N,M) = ZERO
582 A_M(IJK,S,M) = ZERO
583 A_M(IJK,T,M) = ZERO
584 A_M(IJK,B,M) = ZERO
585 A_M(IJK,0,M) = -ONE
586 B_M(IJK,M) = ZERO
587 IF (FLUID_AT(NORTH_OF(IJK))) THEN
588 A_M(IJK,N,M) = -ONE
589 ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
590 A_M(IJK,S,M) = -ONE
591 ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
592 A_M(IJK,T,M) = -ONE
593 ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
594 A_M(IJK,B,M) = -ONE
595 ENDIF
596 ENDDO
597 ENDDO
598 ENDDO
599
600 ELSEIF (BC_TYPE(L) == 'FREE_SLIP_WALL' .AND. .NOT. K_Epsilon) THEN
601 I1 = BC_I_W(L)
602 I2 = BC_I_E(L)
603 J1 = BC_J_S(L)
604 J2 = BC_J_N(L)
605 K1 = BC_K_B(L)
606 K2 = BC_K_T(L)
607 DO K = K1, K2
608 DO J = J1, J2
609 DO I = I1, I2
610 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
611 IF (DEAD_CELL_AT(I,J,K)) CYCLE
612 = FUNIJK(I,J,K)
613 IF (.NOT.WALL_AT(IJK)) CYCLE
614 (IJK,E,M) = ZERO
615 A_M(IJK,W,M) = ZERO
616 A_M(IJK,N,M) = ZERO
617 A_M(IJK,S,M) = ZERO
618 A_M(IJK,T,M) = ZERO
619 A_M(IJK,B,M) = ZERO
620 A_M(IJK,0,M) = -ONE
621 B_M(IJK,M) = ZERO
622 IF (FLUID_AT(NORTH_OF(IJK))) THEN
623 A_M(IJK,N,M) = ONE
624 ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
625 A_M(IJK,S,M) = ONE
626 ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
627 A_M(IJK,T,M) = ONE
628 ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
629 A_M(IJK,B,M) = ONE
630 ENDIF
631 ENDDO
632 ENDDO
633 ENDDO
634
635 ELSEIF (BC_TYPE(L) == 'PAR_SLIP_WALL' .AND. .NOT. K_Epsilon) THEN
636 I1 = BC_I_W(L)
637 I2 = BC_I_E(L)
638 J1 = BC_J_S(L)
639 J2 = BC_J_N(L)
640 K1 = BC_K_B(L)
641 K2 = BC_K_T(L)
642 DO K = K1, K2
643 DO J = J1, J2
644 DO I = I1, I2
645 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
646 IF (DEAD_CELL_AT(I,J,K)) CYCLE
647 = FUNIJK(I,J,K)
648 IF (.NOT.WALL_AT(IJK)) CYCLE
649 = JM1(J)
650 KM = KM1(K)
651 A_M(IJK,E,M) = ZERO
652 A_M(IJK,W,M) = ZERO
653 A_M(IJK,N,M) = ZERO
654 A_M(IJK,S,M) = ZERO
655 A_M(IJK,T,M) = ZERO
656 A_M(IJK,B,M) = ZERO
657 A_M(IJK,0,M) = -ONE
658 B_M(IJK,M) = ZERO
659 IF (FLUID_AT(NORTH_OF(IJK))) THEN
660 IF (BC_HW_G(L) == UNDEFINED) THEN
661 A_M(IJK,N,M) = -HALF
662 A_M(IJK,0,M) = -HALF
663 B_M(IJK,M) = -BC_UW_G(L)
664 ELSE
665 A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODY_N(J))
666 A_M(IJK,N,M) = -(HALF*BC_HW_G(L)-ODY_N(J))
667 B_M(IJK,M) = -BC_HW_G(L)*BC_UW_G(L)
668 ENDIF
669 ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
670 IF (BC_HW_G(L) == UNDEFINED) THEN
671 A_M(IJK,S,M) = -HALF
672 A_M(IJK,0,M) = -HALF
673 B_M(IJK,M) = -BC_UW_G(L)
674 ELSE
675 A_M(IJK,S,M) = -(HALF*BC_HW_G(L)-ODY_N(JM))
676 A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODY_N(JM))
677 B_M(IJK,M) = -BC_HW_G(L)*BC_UW_G(L)
678 ENDIF
679 ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
680 IF (BC_HW_G(L) == UNDEFINED) THEN
681 A_M(IJK,T,M) = -HALF
682 A_M(IJK,0,M) = -HALF
683 B_M(IJK,M) = -BC_UW_G(L)
684 ELSE
685 A_M(IJK,0,M)=-(HALF*BC_HW_G(L)+ODZ_T(K)*OX_E(I))
686 A_M(IJK,T,M)=-(HALF*BC_HW_G(L)-ODZ_T(K)*OX_E(I))
687 B_M(IJK,M) = -BC_HW_G(L)*BC_UW_G(L)
688 ENDIF
689 ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
690 IF (BC_HW_G(L) == UNDEFINED) THEN
691 A_M(IJK,B,M) = -HALF
692 A_M(IJK,0,M) = -HALF
693 B_M(IJK,M) = -BC_UW_G(L)
694 ELSE
695 A_M(IJK,B,M) = -(HALF*BC_HW_G(L)-ODZ_T(KM)*OX_E(I&
696 ))
697 A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODZ_T(KM)*OX_E(I&
698 ))
699 B_M(IJK,M) = -BC_HW_G(L)*BC_UW_G(L)
700 ENDIF
701 ENDIF
702 ENDDO
703 ENDDO
704 ENDDO
705
706
707 ELSEIF (BC_TYPE(L) == 'PAR_SLIP_WALL' .OR. &
708 BC_TYPE(L) == 'NO_SLIP_WALL' .OR. &
709 BC_TYPE(L) == 'FREE_SLIP_WALL' .AND. &
710 K_Epsilon )THEN
711 I1 = BC_I_W(L)
712 I2 = BC_I_E(L)
713 J1 = BC_J_S(L)
714 J2 = BC_J_N(L)
715 K1 = BC_K_B(L)
716 K2 = BC_K_T(L)
717 DO K = K1, K2
718 DO J = J1, J2
719 DO I = I1, I2
720 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
721 IF (DEAD_CELL_AT(I,J,K)) CYCLE
722 = FUNIJK(I,J,K)
723 IF (.NOT.WALL_AT(IJK)) CYCLE
724 = JM1(J)
725 KM = KM1(K)
726 A_M(IJK,E,M) = ZERO
727 A_M(IJK,W,M) = ZERO
728 A_M(IJK,N,M) = ZERO
729 A_M(IJK,S,M) = ZERO
730 A_M(IJK,T,M) = ZERO
731 A_M(IJK,B,M) = ZERO
732 A_M(IJK,0,M) = -ONE
733 B_M(IJK,M) = ZERO
734 IF (FLUID_AT(NORTH_OF(IJK))) THEN
735 CALL Wall_Function(IJK,NORTH_OF(IJK),ODY_N(J),W_F_Slip)
736 A_M(IJK,N,M) = W_F_Slip
737 A_M(IJK,0,M) = -ONE
738 IF (BC_TYPE(L) == 'PAR_SLIP_WALL') B_M(IJK,M) = -BC_UW_G(L)
739 ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
740 CALL Wall_Function(IJK,SOUTH_OF(IJK),ODY_N(JM),W_F_Slip)
741 A_M(IJK,S,M) = W_F_Slip
742 A_M(IJK,0,M) = -ONE
743 IF (BC_TYPE(L) == 'PAR_SLIP_WALL') B_M(IJK,M) = -BC_UW_G(L)
744 ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
745 CALL Wall_Function(IJK,TOP_OF(IJK),ODZ_T(K)*OX_E(I),W_F_Slip)
746 A_M(IJK,T,M) = W_F_Slip
747 A_M(IJK,0,M) = -ONE
748 IF (BC_TYPE(L) == 'PAR_SLIP_WALL') B_M(IJK,M) = -BC_UW_G(L)
749 ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
750 CALL Wall_Function(IJK,BOTTOM_OF(IJK),ODZ_T(KM)*OX_E(I),W_F_Slip)
751 A_M(IJK,B,M) = W_F_Slip
752 A_M(IJK,0,M) = -ONE
753 IF (BC_TYPE(L) == 'PAR_SLIP_WALL') B_M(IJK,M) = -BC_UW_G(L)
754 ENDIF
755 ENDDO
756 ENDDO
757 ENDDO
758
759
760
761
762
763 ELSEIF (BC_TYPE(L)=='P_INFLOW' .OR. BC_TYPE(L)=='P_OUTFLOW') THEN
764 IF (BC_PLANE(L) == 'W') THEN
765
766
767
768 = BC_I_W(L)
769 I2 = BC_I_E(L)
770 J1 = BC_J_S(L)
771 J2 = BC_J_N(L)
772 K1 = BC_K_B(L)
773 K2 = BC_K_T(L)
774 DO K = K1, K2
775 DO J = J1, J2
776 DO I = I1, I2
777 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
778 IF (DEAD_CELL_AT(I,J,K)) CYCLE
779 = FUNIJK(I,J,K)
780 A_M(IJK,E,M) = ZERO
781 A_M(IJK,W,M) = ONE
782 A_M(IJK,N,M) = ZERO
783 A_M(IJK,S,M) = ZERO
784 A_M(IJK,T,M) = ZERO
785 A_M(IJK,B,M) = ZERO
786 A_M(IJK,0,M) = -ONE
787 B_M(IJK,M) = ZERO
788 ENDDO
789 ENDDO
790 ENDDO
791 ENDIF
792
793
794
795
796
797 ELSEIF (BC_TYPE(L) == 'OUTFLOW') THEN
798 IF (BC_PLANE(L) == 'W') THEN
799 I1 = BC_I_W(L)
800 I2 = BC_I_E(L)
801 J1 = BC_J_S(L)
802 J2 = BC_J_N(L)
803 K1 = BC_K_B(L)
804 K2 = BC_K_T(L)
805 DO K = K1, K2
806 DO J = J1, J2
807 DO I = I1, I2
808 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
809 IF (DEAD_CELL_AT(I,J,K)) CYCLE
810 = FUNIJK(I,J,K)
811 A_M(IJK,E,M) = ZERO
812 A_M(IJK,W,M) = ONE
813 A_M(IJK,N,M) = ZERO
814 A_M(IJK,S,M) = ZERO
815 A_M(IJK,T,M) = ZERO
816 A_M(IJK,B,M) = ZERO
817 A_M(IJK,0,M) = -ONE
818 B_M(IJK,M) = ZERO
819 IM = IM1(I)
820 IMJK = IM_OF(IJK)
821 A_M(IMJK,E,M) = ZERO
822 A_M(IMJK,W,M) = X_E(IM)/X_E(IM1(IM))
823 A_M(IMJK,N,M) = ZERO
824 A_M(IMJK,S,M) = ZERO
825 A_M(IMJK,T,M) = ZERO
826 A_M(IMJK,B,M) = ZERO
827 A_M(IMJK,0,M) = -ONE
828 B_M(IMJK,M) = ZERO
829 ENDDO
830 ENDDO
831 ENDDO
832 ELSEIF (BC_PLANE(L) == 'E') THEN
833 I1 = BC_I_W(L)
834 I2 = BC_I_E(L)
835 J1 = BC_J_S(L)
836 J2 = BC_J_N(L)
837 K1 = BC_K_B(L)
838 K2 = BC_K_T(L)
839 DO K = K1, K2
840 DO J = J1, J2
841 DO I = I1, I2
842 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
843 IF (DEAD_CELL_AT(I,J,K)) CYCLE
844 = FUNIJK(I,J,K)
845 IP = IP1(I)
846 IPJK = IP_OF(IJK)
847 A_M(IPJK,E,M) = X_E(IP)/X_E(I)
848 A_M(IPJK,W,M) = ZERO
849 A_M(IPJK,N,M) = ZERO
850 A_M(IPJK,S,M) = ZERO
851 A_M(IPJK,T,M) = ZERO
852 A_M(IPJK,B,M) = ZERO
853 A_M(IPJK,0,M) = -ONE
854 B_M(IPJK,M) = ZERO
855 ENDDO
856 ENDDO
857 ENDDO
858 ENDIF
859
860
861
862
863
864
865
866 ELSE
867 I1 = BC_I_W(L)
868 I2 = BC_I_E(L)
869 J1 = BC_J_S(L)
870 J2 = BC_J_N(L)
871 K1 = BC_K_B(L)
872 K2 = BC_K_T(L)
873 DO K = K1, K2
874 DO J = J1, J2
875 DO I = I1, I2
876 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
877 IF (DEAD_CELL_AT(I,J,K)) CYCLE
878 = FUNIJK(I,J,K)
879
880 (IJK,E,M) = ZERO
881 A_M(IJK,W,M) = ZERO
882 A_M(IJK,N,M) = ZERO
883 A_M(IJK,S,M) = ZERO
884 A_M(IJK,T,M) = ZERO
885 A_M(IJK,B,M) = ZERO
886 A_M(IJK,0,M) = -ONE
887 B_M(IJK,M) = -U_G(IJK)
888 IF (BC_PLANE(L) == 'W') THEN
889
890
891
892 = WEST_OF(IJK)
893 A_M(IJKW,E,M) = ZERO
894 A_M(IJKW,W,M) = ZERO
895 A_M(IJKW,N,M) = ZERO
896 A_M(IJKW,S,M) = ZERO
897 A_M(IJKW,T,M) = ZERO
898 A_M(IJKW,B,M) = ZERO
899 A_M(IJKW,0,M) = -ONE
900 B_M(IJKW,M) = -U_G(IJKW)
901 ENDIF
902 ENDDO
903 ENDDO
904 ENDDO
905 ENDIF
906
907
908
909
910
911
912 ENDIF
913 ENDDO
914
915 RETURN
916 END SUBROUTINE SOURCE_U_G_BC
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937 SUBROUTINE Wall_Function(IJK1, IJK2, ODX_WF, W_F_Slip)
938
939
940
941
942 USE param
943 USE param1
944 USE physprop
945 USE fldvar
946 USE visc_g
947 USE geometry
948 USE indices
949 USE bc
950 USE compar
951 USE turb
952 USE mpi_utility
953 IMPLICIT NONE
954
955
956
957
958 INTEGER :: IJK1, IJK2
959
960 DOUBLE PRECISION ODX_WF, W_F_Slip
961
962
963
964
965 DOUBLE PRECISION, PARAMETER :: C_mu = 0.09D0
966
967 DOUBLE PRECISION, PARAMETER :: Kappa = 0.42D0
968
969
970
971
972
973
974 IF(DABS(ODX_WF)>1.0D-5) THEN
975
976
977
978 = (ONE - ONE/ODX_WF* RO_g(IJK2)*C_mu**0.25 * &
979 SQRT(K_Turb_G(IJK2))/MU_gT(IJK2) * &
980 Kappa/LOG(9.81D+0/(ODX_WF*2.D+0)*RO_g(IJK2)*C_mu**0.25*&
981 SQRT(K_Turb_G(IJK2))/MU_g(IJK2)))
982 ELSE
983
984 = ONE
985 ENDIF
986
987
988 RETURN
989 END SUBROUTINE Wall_Function
990
991
992
993
994
995
996
997
998
999
1000
1001
1002 SUBROUTINE POINT_SOURCE_U_G(A_M, B_M)
1003
1004 use compar
1005 use constant
1006 use geometry
1007 use indices
1008 use physprop
1009 use ps
1010 use run
1011 use functions
1012
1013 IMPLICIT NONE
1014
1015
1016
1017
1018 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
1019
1020 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
1021
1022
1023
1024
1025 INTEGER :: IJK, I, J, K
1026 INTEGER :: PSV, M
1027 INTEGER :: lIE, lIW
1028
1029 DOUBLE PRECISION :: pSource
1030
1031
1032
1033 = 0
1034
1035
1036
1037 PS_LP: do PSV = 1, DIMENSION_PS
1038 if(.NOT.PS_DEFINED(PSV)) cycle PS_LP
1039 if(abs(PS_U_g(PSV)) < small_number) cycle PS_LP
1040
1041 if(PS_U_g(PSV) < 0.0d0) then
1042 lIW = PS_I_W(PSV) - 1
1043 lIE = PS_I_E(PSV) - 1
1044 else
1045 lIW = PS_I_W(PSV)
1046 lIE = PS_I_E(PSV)
1047 endif
1048
1049 do k = PS_K_B(PSV), PS_K_T(PSV)
1050 do j = PS_J_S(PSV), PS_J_N(PSV)
1051 do i = lIW, lIE
1052
1053 if(.NOT.IS_ON_myPE_plus2layers(I,J,K)) cycle
1054 IF (DEAD_CELL_AT(I,J,K)) CYCLE
1055
1056 = funijk(i,j,k)
1057 if(.NOT. fluid_at(ijk)) cycle
1058
1059 pSource = PS_MASSFLOW_G(PSV) * (VOL(IJK)/PS_VOLUME(PSV))
1060
1061 B_M(IJK,M) = B_M(IJK,M) - pSource * &
1062 PS_U_g(PSV) * PS_VEL_MAG_g(PSV)
1063
1064 enddo
1065 enddo
1066 enddo
1067 enddo PS_LP
1068
1069 RETURN
1070 END SUBROUTINE POINT_SOURCE_U_G
1071