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