File: N:\mfix\model\cartesian_grid\get_alpha.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 SUBROUTINE GET_3D_ALPHA_U_CUT_CELL
16
17 USE bc
18 USE compar, ONLY: iend1, jend1, kend1, IJKSTART3, IJKEND3, mype, pe_io
19 USE cutcell
20 USE exit, only: mfix_exit
21 USE functions, ONLY: FUNIJK
22 USE geometry, ONLY: DO_K, NO_K, ayz, ayz_u, flag_e
23 USE indices, ONLY: I_OF, J_OF, K_OF
24 USE param1, ONLY: half, one, undefined, zero
25
26 IMPLICIT NONE
27 DOUBLE PRECISION:: Xe,Ye,Ze,Xn,Yn,Zn,Xt,Yt,Zt
28 INTEGER :: I,J,K,IM,JM,KM,IP,JP,KP,IJK,IMJK,IJMK,IPJK,IJPK,IJKM,IJKP
29 DOUBLE PRECISION :: Xi
30 DOUBLE PRECISION :: Sx,Sy,Sz
31 DOUBLE PRECISION :: Nx,Ny,Nz
32 DOUBLE PRECISION :: DELH_ec , DELH_e
33 DOUBLE PRECISION :: DELH_nc , DELH_n
34 DOUBLE PRECISION :: DELH_tc , DELH_t
35 LOGICAL :: V_NODE_AT_NE,V_NODE_AT_NW
36 LOGICAL :: W_NODE_AT_TE,W_NODE_AT_TW
37 INTEGER :: BCV
38 INTEGER ::BCT
39
40 IF(MyPE == PE_IO) THEN
41 WRITE(*,10)'COMPUTING INTERPOLATION FACTORS IN U-MOMENTUM CELLS...'
42 ENDIF
43 10 FORMAT(1X,A)
44
45 DELH_U = UNDEFINED
46
47 Theta_Ue = HALF
48 Theta_Ue_bar = HALF
49 Theta_U_ne = HALF
50 Theta_U_nw = HALF
51 Theta_U_te = HALF
52 Theta_U_tw = HALF
53 ALPHA_Ue_c = ONE
54 NOC_U_E = ZERO
55 Theta_Un = HALF
56 Theta_Un_bar = HALF
57 ALPHA_Un_c = ONE
58 NOC_U_N = ZERO
59 Theta_Ut = HALF
60 Theta_Ut_bar = HALF
61 ALPHA_Ut_c = ONE
62 NOC_U_T = ZERO
63
64 A_UPG_E = AYZ_U
65 A_UPG_W = AYZ_U
66
67
68 DO IJK = IJKSTART3, IJKEND3
69 IF(INTERIOR_CELL_AT(IJK)) THEN
70
71 I = I_OF(IJK)
72 J = J_OF(IJK)
73 K = K_OF(IJK)
74
75 IM = I - 1
76 JM = J - 1
77 KM = K - 1
78
79 IP = I + 1
80 JP = J + 1
81 KP = K + 1
82
83 IMJK = FUNIJK(IM,J,K)
84 IPJK = FUNIJK(IP,J,K)
85 IJMK = FUNIJK(I,JM,K)
86 IJPK = FUNIJK(I,JP,K)
87 IJKM = FUNIJK(I,J,KM)
88 IJKP = FUNIJK(I,J,KP)
89
90 CALL GET_CELL_NODE_COORDINATES(IJK,'U_MOMENTUM')
91
92
93
94
95
96 (IJK) = DELX_Ue(IJK) / (DELX_Ue(IJK) + DELX_Uw(IPJK))
97 Theta_Ue_bar(IJK) = ONE - Theta_Ue(IJK)
98
99
100
101
102
103 = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.(.NOT.WALL_V_AT(IJK)))
104 V_NODE_AT_NE = ((.NOT.BLOCKED_V_CELL_AT(IPJK)).AND.(.NOT.WALL_V_AT(IPJK)))
105
106 IF(V_NODE_AT_NW.AND.V_NODE_AT_NE) THEN
107 Theta_U_ne(IJK) = (X_U_nc(IJK) - X_V(IJK) ) / (X_V(IPJK) - X_V(IJK))
108 Theta_U_nw(IJK) = ONE - Theta_U_ne(IJK)
109 ELSE IF (V_NODE_AT_NE.AND.(.NOT.V_NODE_AT_NW)) THEN
110 IF(NO_K) THEN
111 Xi = Xn_U_int(IJK)
112 ELSE
113 Xi = HALF * (Xn_U_int(IJK) + Xn_U_int(IJKM))
114 ENDIF
115 Theta_U_ne(IJK) = (X_U_nc(IJK) - Xi) / (X_V(IPJK) - Xi)
116 Theta_U_nw(IJK) = ONE - Theta_U_ne(IJK)
117 ELSE IF ((.NOT.V_NODE_AT_NE).AND.V_NODE_AT_NW) THEN
118 IF(NO_K) THEN
119 Xi = Xn_U_int(IJK)
120 ELSE
121 Xi = HALF * (Xn_U_int(IJK) + Xn_U_int(IJKM))
122 ENDIF
123 Theta_U_ne(IJK) = (X_U_nc(IJK) - X_V(IJK) ) / (Xi - X_V(IJK))
124 Theta_U_nw(IJK) = ONE - Theta_U_ne(IJK)
125 ELSE
126 Theta_U_ne(IJK) = ZERO
127 Theta_U_nw(IJK) = ZERO
128 ENDIF
129
130
131 IF ((Theta_U_ne(IJK)>=ONE).OR.(Theta_U_ne(IJK)<=ZERO)) THEN
132 Theta_U_ne(IJK) = HALF
133 Theta_U_nw(IJK) = HALF
134 ENDIF
135
136
137
138
139
140 IF(DO_K) THEN
141
142 W_NODE_AT_TW = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.(.NOT.WALL_W_AT(IJK)))
143 W_NODE_AT_TE = ((.NOT.BLOCKED_W_CELL_AT(IPJK)).AND.(.NOT.WALL_W_AT(IPJK)))
144
145 IF(W_NODE_AT_TW.AND.W_NODE_AT_TE) THEN
146 Theta_U_te(IJK) = (X_U_tc(IJK) - X_W(IJK) ) / (X_W(IPJK) - X_W(IJK))
147 Theta_U_tw(IJK) = ONE - Theta_U_te(IJK)
148
149 ELSE IF (W_NODE_AT_TE.AND.(.NOT.W_NODE_AT_TW)) THEN
150 Xi = HALF * (Xn_U_int(IJK) + Xn_U_int(IJMK))
151 Theta_U_te(IJK) = (X_U_tc(IJK) - Xi) / (X_W(IPJK) - Xi)
152 Theta_U_tw(IJK) = ONE - Theta_U_te(IJK)
153 ELSE IF ((.NOT.W_NODE_AT_TE).AND.W_NODE_AT_TW) THEN
154 Xi = HALF * (Xn_U_int(IJK) + Xn_U_int(IJMK))
155 Theta_U_te(IJK) = (X_U_tc(IJK) - X_W(IJK) ) / (Xi - X_W(IJK))
156 Theta_U_tw(IJK) = ONE - Theta_U_te(IJK)
157 ELSE
158 Theta_U_te(IJK) = ZERO
159 Theta_U_tw(IJK) = ZERO
160 ENDIF
161
162 IF ((Theta_U_te(IJK)>=ONE).OR.(Theta_U_te(IJK)<=ZERO)) THEN
163 Theta_U_te(IJK) = HALF
164 Theta_U_tw(IJK) = HALF
165 ENDIF
166 ENDIF
167
168 BCV = BC_U_ID(IJK)
169
170 IF(BCV>0) THEN
171 BCT = BC_TYPE_ENUM(BCV)
172 ELSE
173 BCT = BLANK
174 ENDIF
175
176 IF(BCT ==CG_NSW.OR.BCT ==CG_PSW) THEN
177
178 CALL GET_DEL_H(IJK,'U_MOMENTUM',X_U(IJK),Y_U(IJK),Z_U(IJK),DELH_U(IJK),Nx,Ny,Nz)
179
180 IF(DELH_U(IJK)<ZERO.AND.(.NOT.WALL_U_AT(IJK))) THEN
181 WRITE(*,*) 'NEGATIVE DELH-U AT XYZ=', X_U(IJK),Y_U(IJK),Z_U(IJK)
182 WRITE(*,*) 'DELH_U=', DELH_U(IJK)
183 WRITE(*,*) 'AYZ=', AYZ(IJK)
184 WRITE(*,*) 'MFIX WILL EXIT NOW.'
185 CALL MFIX_EXIT(MYPE)
186 ENDIF
187
188
189
190
191
192
193
194
195 = X_NODE(8)
196 Ye = Theta_Ue_bar(IJK) * Y_U(IJK) + Theta_Ue(IJK) * Y_U(IPJK)
197 Ze = Theta_Ue_bar(IJK) * Z_U(IJK) + Theta_Ue(IJK) * Z_U(IPJK)
198
199 CALL GET_DEL_H(IJK,'U_MOMENTUM',X_U_ec(IJK),Y_U_ec(IJK),Z_U_ec(IJK),DELH_ec,Nx,Ny,Nz)
200 CALL GET_DEL_H(IJK,'U_MOMENTUM',Xe,Ye,Ze,DELH_e,Nx,Ny,Nz)
201
202 IF((DELH_ec == UNDEFINED).OR.(DELH_e == UNDEFINED)) THEN
203 ALPHA_Ue_c(IJK) = ONE
204 ELSE
205 ALPHA_Ue_c(IJK) = DMIN1(ALPHA_MAX , DELH_ec / DELH_e )
206 ENDIF
207
208
209 IF(BLOCKED_U_CELL_AT(IPJK)) ALPHA_Ue_c(IJK) = ZERO
210 IF(WALL_U_AT(IJK).AND.WALL_U_AT(IPJK)) ALPHA_Ue_c(IJK) = ZERO
211 IF(I == IEND1) ALPHA_Ue_c(IJK) = ONE
212
213 IF(ALPHA_Ue_c(IJK)<ZERO) THEN
214 WRITE(*,*) 'NEGATIVE ALPHA_Ue_c at IJK=',IJK
215 WRITE(*,*) 'MFIX WILL EXIT NOW.'
216 CALL MFIX_EXIT(MYPE)
217 ENDIF
218
219
220
221
222
223 = X_U(IPJK) - X_U(IJK)
224 Sy = Y_U(IPJK) - Y_U(IJK)
225 Sz = Z_U(IPJK) - Z_U(IJK)
226
227 NOC_U_E(IJK) = (Sy * Ny + Sz * Nz)/(Sx * DELH_e)
228
229 IF(BLOCKED_U_CELL_AT(IPJK)) NOC_U_E(IJK) = ZERO
230 IF(WALL_U_AT(IJK).AND.WALL_U_AT(IPJK)) NOC_U_E(IJK) = ZERO
231 IF(I == IEND1) NOC_U_E(IJK) = ZERO
232
233
234
235
236
237 (IJK) = DELY_Un(IJK) / (DELY_Un(IJK) + DELY_Us(IJPK))
238 Theta_Un_bar(IJK) = ONE - Theta_Un(IJK)
239
240 IF ((Theta_Un(IJK)>=ONE).OR.(Theta_Un(IJK)<=ZERO)) THEN
241 Theta_Un(IJK) = HALF
242 Theta_Un_bar(IJK) = HALF
243 ENDIF
244
245
246
247
248
249
250 = Theta_Un_bar(IJK) * X_U(IJK) + Theta_Un(IJK) * X_U(IJPK)
251 Yn = Y_NODE(8)
252 Zn = Theta_Un_bar(IJK) * Z_U(IJK) + Theta_Un(IJK) * Z_U(IJPK)
253
254 CALL GET_DEL_H(IJK,'U_MOMENTUM',X_U_nc(IJK),Y_U_nc(IJK),Z_U_nc(IJK),DELH_nc,Nx,Ny,Nz)
255 CALL GET_DEL_H(IJK,'U_MOMENTUM',Xn,Yn,Zn,DELH_n,Nx,Ny,Nz)
256
257 IF((DELH_nc == UNDEFINED).OR.(DELH_n == UNDEFINED)) THEN
258 ALPHA_Un_c(IJK) = ONE
259 ELSE
260 ALPHA_Un_c(IJK) = DMIN1(ALPHA_MAX , DELH_nc / DELH_n )
261 ENDIF
262
263
264 IF(BLOCKED_U_CELL_AT(IJPK)) ALPHA_Un_c(IJK) = ZERO
265 IF(WALL_U_AT(IJK).AND.WALL_U_AT(IJPK)) ALPHA_Un_c(IJK) = ZERO
266 IF(J == JEND1) ALPHA_Un_c(IJK) = ONE
267
268 IF(ALPHA_Un_c(IJK)<ZERO) ALPHA_Un_c(IJK) = ONE
269
270
271
272
273
274 = X_U(IJPK) - X_U(IJK)
275 Sy = Y_U(IJPK) - Y_U(IJK)
276 Sz = Z_U(IJPK) - Z_U(IJK)
277
278 NOC_U_N(IJK) = (Sx * Nx + Sz * Nz)/(Sy * DELH_n)
279
280 IF(BLOCKED_U_CELL_AT(IJPK)) NOC_U_N(IJK) = ZERO
281 IF(WALL_U_AT(IJK).AND.WALL_U_AT(IJPK)) NOC_U_N(IJK) = ZERO
282 IF(J == JEND1) NOC_U_N(IJK) = ZERO
283
284
285 IF(DO_K) THEN
286
287
288
289
290 (IJK) = DELZ_Ut(IJK) / (DELZ_Ut(IJK) + DELZ_Ub(IJKP))
291 Theta_Ut_bar(IJK) = ONE - Theta_Ut(IJK)
292
293 IF ((Theta_Ut(IJK)>=ONE).OR.(Theta_Ut(IJK)<=ZERO)) THEN
294 Theta_Ut(IJK) = HALF
295 Theta_Ut_bar(IJK) = HALF
296 ENDIF
297
298
299
300
301
302 = Theta_Ut_bar(IJK) * X_U(IJK) + Theta_Ut(IJK) * X_U(IJKP)
303 Yt = Theta_Ut_bar(IJK) * Y_U(IJK) + Theta_Ut(IJK) * Y_U(IJKP)
304 Zt = Z_NODE(8)
305
306 CALL GET_DEL_H(IJK,'U_MOMENTUM',X_U_tc(IJK),Y_U_tc(IJK),Z_U_tc(IJK),DELH_tc,Nx,Ny,Nz)
307 CALL GET_DEL_H(IJK,'U_MOMENTUM',Xt,Yt,Zt,DELH_t,Nx,Ny,Nz)
308
309 IF((DELH_tc == UNDEFINED).OR.(DELH_t == UNDEFINED)) THEN
310 ALPHA_Ut_c(IJK) = ONE
311 ELSE
312 ALPHA_Ut_c(IJK) = DMIN1(ALPHA_MAX , DELH_tc / DELH_t )
313 ENDIF
314
315 IF(BLOCKED_U_CELL_AT(IJKP)) ALPHA_Ut_c(IJK) = ZERO
316 IF(WALL_U_AT(IJK).AND.WALL_U_AT(IJKP)) ALPHA_Ut_c(IJK) = ZERO
317 IF(K == KEND1) ALPHA_Ut_c(IJK) = ONE
318
319 IF(ALPHA_Ut_c(IJK)<ZERO) ALPHA_Ut_c(IJK) = ONE
320
321
322
323
324
325 = X_U(IJKP) - X_U(IJK)
326 Sy = Y_U(IJKP) - Y_U(IJK)
327 Sz = Z_U(IJKP) - Z_U(IJK)
328
329 NOC_U_T(IJK) = (Sx * Nx + Sy * Ny)/(Sz * DELH_t)
330
331
332 IF(BLOCKED_U_CELL_AT(IJKP)) NOC_U_T(IJK) = ZERO
333 IF(WALL_U_AT(IJK).AND.WALL_U_AT(IJKP)) NOC_U_T(IJK) = ZERO
334 IF(K == KEND1) NOC_U_T(IJK) = ZERO
335
336 ENDIF
337
338 ENDIF
339
340
341
342
343
344
345
346 IF ( CUT_U_CELL_AT(IJK) ) THEN
347
348 SELECT CASE (PG_OPTION)
349
350 CASE(0)
351
352
353 CASE(1)
354
355 A_UPG_E(IJK) = dmax1(AYZ_U(IJK),AYZ_U(IMJK))
356 A_UPG_W(IJK) = A_UPG_E(IJK)
357
358 CASE(2)
359
360 A_UPG_E(IJK) = AYZ_U(IJK)
361 A_UPG_W(IJK) = AYZ_U(IMJK)
362
363 CASE DEFAULT
364
365 WRITE(*,*)'INVALID PRESSURE GRADIENT OPTION:',PG_OPTION
366 WRITE(*,*)'PG_OPTION SHOULD BE SET EQUAL TO 0, 1 OR 2.'
367 WRITE(*,*)'PLEASE VERIFY MFIX.DAT FILE AND TRY AGAIN.'
368 WRITE(*,*) 'MFIX WILL EXIT NOW.'
369 CALL MFIX_EXIT(myPE)
370
371 END SELECT
372
373
374 IF (BLOCKED_CELL_AT(IJK).OR.BLOCKED_CELL_AT(IPJK)) THEN
375 IF(.NOT.WALL_U_AT(IJK)) THEN
376 WALL_U_AT(IJK) = .TRUE.
377 FLAG_E(IJK) = 0
378 IF(PRINT_WARNINGS) THEN
379 WRITE(*,*) 'WARNING: ONLY ONE PRESSURE NODE DETECTED IN U-MOMENTUM CELL IJK =',IJK
380 WRITE(*,*) 'RESETTING U-MOMENTUM CELL AS U_WALL CELL.'
381 ENDIF
382
383
384 ENDIF
385 ENDIF
386
387 ELSE
388
389 (IJK) = AYZ_U(IJK)
390 A_UPG_W(IJK) = AYZ_U(IJK)
391
392 ENDIF
393
394 ENDIF
395 END DO
396
397 IF(PG_OPTION==0) THEN
398 A_UPG_E = AYZ
399 A_UPG_W = AYZ
400 ENDIF
401
402 RETURN
403
404 END SUBROUTINE GET_3D_ALPHA_U_CUT_CELL
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420 SUBROUTINE GET_3D_ALPHA_V_CUT_CELL
421
422 USE bc
423 USE compar, ONLY: iend1, jend1, kend1, IJKSTART3, IJKEND3, mype, pe_io
424 USE cutcell
425 USE functions, ONLY: FUNIJK
426 USE exit, only: mfix_exit
427 USE geometry, ONLY: DO_K, NO_K, axz, axz_v, flag_n
428 USE indices, ONLY: I_OF, J_OF, K_OF
429 USE param1, ONLY: half, one, undefined, zero
430
431 IMPLICIT NONE
432 DOUBLE PRECISION:: Xe,Ye,Ze,Xn,Yn,Zn,Xt,Yt,Zt
433 INTEGER :: I,J,K,IM,JM,KM,IP,JP,KP,IJK,IMJK,IJMK,IPJK,IJPK,IJKM,IJKP
434 DOUBLE PRECISION :: Yi
435 DOUBLE PRECISION :: Sx,Sy,Sz
436 DOUBLE PRECISION :: Nx,Ny,Nz
437 DOUBLE PRECISION :: DELH_ec , DELH_e
438 DOUBLE PRECISION :: DELH_nc , DELH_n
439 DOUBLE PRECISION :: DELH_tc , DELH_t
440 LOGICAL :: U_NODE_AT_NE, U_NODE_AT_SE
441 LOGICAL :: W_NODE_AT_NT, W_NODE_AT_ST
442 INTEGER :: BCV
443 INTEGER ::BCT
444
445 IF(MyPE == PE_IO) THEN
446 WRITE(*,10)'COMPUTING INTERPOLATION FACTORS IN V-MOMENTUM CELLS...'
447 ENDIF
448 10 FORMAT(1X,A)
449
450 DELH_V = UNDEFINED
451
452 Theta_V_ne = HALF
453 Theta_V_se = HALF
454 Theta_Vn = HALF
455 Theta_Vn_bar = HALF
456 Theta_V_nt = HALF
457 Theta_V_st = HALF
458 Theta_Ve = HALF
459 Theta_Ve_bar = HALF
460 ALPHA_Ve_c = ONE
461 NOC_V_E = ZERO
462 ALPHA_Vn_c = ONE
463 NOC_V_N = ZERO
464 Theta_Vt = HALF
465 Theta_Vt_bar = HALF
466 ALPHA_Vt_c = ONE
467 NOC_V_T = ZERO
468
469 A_VPG_N = AXZ_V
470 A_VPG_S = AXZ_V
471
472
473 DO IJK = IJKSTART3, IJKEND3
474 IF(INTERIOR_CELL_AT(IJK)) THEN
475
476 I = I_OF(IJK)
477 J = J_OF(IJK)
478 K = K_OF(IJK)
479
480 IM = I - 1
481 JM = J - 1
482 KM = K - 1
483
484 IP = I + 1
485 JP = J + 1
486 KP = K + 1
487
488 IMJK = FUNIJK(IM,J,K)
489 IPJK = FUNIJK(IP,J,K)
490 IJMK = FUNIJK(I,JM,K)
491 IJPK = FUNIJK(I,JP,K)
492 IJKM = FUNIJK(I,J,KM)
493 IJKP = FUNIJK(I,J,KP)
494
495 CALL GET_CELL_NODE_COORDINATES(IJK,'V_MOMENTUM')
496
497
498
499
500
501 = ((.NOT.BLOCKED_U_CELL_AT(IJPK)).AND.(.NOT.WALL_U_AT(IJPK)))
502 U_NODE_AT_SE = ((.NOT.BLOCKED_U_CELL_AT(IJK)).AND.(.NOT.WALL_U_AT(IJK)))
503
504 IF(U_NODE_AT_SE.AND.U_NODE_AT_NE) THEN
505 Theta_V_ne(IJK) = (Y_V_ec(IJK) - Y_U(IJK) ) / (Y_U(IJPK) - Y_U(IJK))
506 Theta_V_se(IJK) = ONE - Theta_V_ne(IJK)
507 ELSE IF (U_NODE_AT_SE.AND.(.NOT.U_NODE_AT_NE)) THEN
508 IF(NO_K) THEN
509 Yi = Ye_V_int(IJK)
510 ELSE
511 Yi = HALF * (Ye_V_int(IJK) + Ye_V_int(IJKM))
512 ENDIF
513 Theta_V_ne(IJK) = (Y_V_ec(IJK) - Y_U(IJK) ) / (Yi - Y_U(IJK))
514 Theta_V_se(IJK) = ONE - Theta_V_ne(IJK)
515 ELSE IF ((.NOT.U_NODE_AT_SE).AND.U_NODE_AT_NE) THEN
516 IF(NO_K) THEN
517 Yi = Ye_V_int(IJK)
518 ELSE
519 Yi = HALF * (Ye_V_int(IJK) + Ye_V_int(IJKM))
520 ENDIF
521 Theta_V_ne(IJK) = (Y_V_ec(IJK) - Yi ) / (Y_U(IJPK) - Yi)
522 Theta_V_se(IJK) = ONE - Theta_V_ne(IJK)
523 ELSE
524 Theta_V_ne(IJK) = ZERO
525 Theta_V_se(IJK) = ZERO
526 ENDIF
527
528 IF ((Theta_V_ne(IJK)>=ONE).OR.(Theta_V_ne(IJK)<=ZERO)) THEN
529 Theta_V_ne(IJK) = HALF
530 Theta_V_se(IJK) = HALF
531 ENDIF
532
533
534
535
536
537 (IJK) = DELY_Vn(IJK) / (DELY_Vn(IJK) + DELY_Vs(IJPK))
538 Theta_Vn_bar(IJK) = ONE - Theta_Vn(IJK)
539
540 IF ((Theta_Vn(IJK)>=ONE).OR.(Theta_Vn(IJK)<=ZERO)) THEN
541 Theta_Vn(IJK) = HALF
542 Theta_Vn_bar(IJK) = HALF
543 ENDIF
544
545 IF(DO_K) THEN
546
547
548
549
550 = ((.NOT.BLOCKED_W_CELL_AT(IJPK)).AND.(.NOT.WALL_W_AT(IJPK)))
551 W_NODE_AT_ST = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.(.NOT.WALL_W_AT(IJK)))
552
553 IF(W_NODE_AT_ST.AND.W_NODE_AT_NT) THEN
554 Theta_V_nt(IJK) = (Y_V_tc(IJK) - Y_W(IJK) ) / (Y_W(IJPK) - Y_W(IJK))
555 Theta_V_st(IJK) = ONE - Theta_V_nt(IJK)
556 ELSE IF (W_NODE_AT_ST.AND.(.NOT.W_NODE_AT_NT)) THEN
557 Yi = HALF * (Ye_V_int(IJK) + Ye_V_int(IMJK))
558 Theta_V_nt(IJK) = (Y_V_tc(IJK) - Y_W(IJK) ) / (Yi - Y_W(IJK))
559 Theta_V_st(IJK) = ONE - Theta_V_nt(IJK)
560 ELSE IF ((.NOT.W_NODE_AT_ST).AND.W_NODE_AT_NT) THEN
561 Yi = HALF * (Ye_V_int(IJK) + Ye_V_int(IMJK))
562 Theta_V_nt(IJK) = (Y_V_tc(IJK) - Yi ) / (Y_W(IJPK) - Yi)
563 Theta_V_st(IJK) = ONE - Theta_V_nt(IJK)
564 ELSE
565 Theta_V_nt(IJK) = ZERO
566 Theta_V_st(IJK) = ZERO
567 ENDIF
568
569 IF ((Theta_V_nt(IJK)>=ONE).OR.(Theta_V_nt(IJK)<=ZERO)) THEN
570 Theta_V_nt(IJK) = HALF
571 Theta_V_st(IJK) = HALF
572 ENDIF
573
574 ENDIF
575
576 BCV = BC_V_ID(IJK)
577 IF(BCV>0) THEN
578 BCT = BC_TYPE_ENUM(BCV)
579 ELSE
580 BCT = BLANK
581 ENDIF
582
583 IF(BCT ==CG_NSW.OR.BCT ==CG_PSW) THEN
584
585 CALL GET_DEL_H(IJK,'V_MOMENTUM',X_V(IJK),Y_V(IJK),Z_V(IJK),DELH_V(IJK),Nx,Ny,Nz)
586
587 IF(DELH_V(IJK)<ZERO.AND.(.NOT.WALL_V_AT(IJK))) THEN
588 WRITE(*,*) 'NEGATIVE DELH-V AT XYZ=', X_V(IJK),Y_V(IJK),Z_V(IJK)
589 WRITE(*,*) 'DELH_V=', DELH_V(IJK)
590 WRITE(*,*) 'AYZ=', AXZ(IJK)
591 WRITE(*,*) 'MFIX WILL EXIT NOW.'
592 CALL MFIX_EXIT(MYPE)
593 ENDIF
594
595
596
597
598
599 (IJK) = DELX_Ve(IJK) / (DELX_Ve(IJK) + DELX_Vw(IPJK))
600 Theta_Ve_bar(IJK) = ONE - Theta_Ve(IJK)
601
602 IF ((Theta_Ve(IJK)>=ONE).OR.(Theta_Ve(IJK)<=ZERO)) THEN
603 Theta_Ve(IJK) = HALF
604 Theta_Ve_bar(IJK) = HALF
605 ENDIF
606
607
608
609
610 = X_NODE(8)
611 Ye = Theta_Ve_bar(IJK) * Y_V(IJK) + Theta_Ve(IJK) * Y_V(IPJK)
612 Ze = Theta_Ve_bar(IJK) * Z_V(IJK) + Theta_Ve(IJK) * Z_V(IPJK)
613
614
615 CALL GET_DEL_H(IJK,'V_MOMENTUM',X_V_ec(IJK),Y_V_ec(IJK),Z_V_ec(IJK),DELH_ec,Nx,Ny,Nz)
616
617 CALL GET_DEL_H(IJK,'V_MOMENTUM',Xe,Ye,Ze,DELH_e,Nx,Ny,Nz)
618
619 IF((DELH_ec == UNDEFINED).OR.(DELH_e == UNDEFINED)) THEN
620 ALPHA_Ve_c(IJK) = ONE
621 ELSE
622 ALPHA_Ve_c(IJK) = DMIN1(ALPHA_MAX , DELH_ec / DELH_e )
623 ENDIF
624
625
626 IF(BLOCKED_V_CELL_AT(IPJK)) ALPHA_Ve_c(IJK) = ZERO
627 IF(WALL_V_AT(IJK).AND.WALL_V_AT(IPJK)) ALPHA_Ve_c(IJK) = ZERO
628 IF(I == IEND1) ALPHA_Ve_c(IJK) = ONE
629
630 IF(ALPHA_Ve_c(IJK)<ZERO) ALPHA_Ve_c(IJK) = ONE
631
632
633
634
635
636 = X_V(IPJK) - X_V(IJK)
637 Sy = Y_V(IPJK) - Y_V(IJK)
638 Sz = Z_V(IPJK) - Z_V(IJK)
639
640 NOC_V_E(IJK) = (Sy * Ny + Sz * Nz)/(Sx * DELH_e)
641
642 IF(BLOCKED_V_CELL_AT(IPJK)) NOC_V_E(IJK) = ZERO
643 IF(WALL_V_AT(IJK).AND.WALL_V_AT(IPJK)) NOC_V_E(IJK) = ZERO
644 IF(I == IEND1) NOC_V_E(IJK) = ZERO
645
646
647
648
649
650
651
652
653 = Theta_Vn_bar(IJK) * X_V(IJK) + Theta_Vn(IJK) * X_V(IJPK)
654 Yn = Y_NODE(8)
655 Zn = Theta_Vn_bar(IJK) * Z_V(IJK) + Theta_Vn(IJK) * Z_V(IJPK)
656
657 CALL GET_DEL_H(IJK,'V_MOMENTUM',X_V_nc(IJK),Y_V_nc(IJK),Z_V_nc(IJK),DELH_nc,Nx,Ny,Nz)
658
659 CALL GET_DEL_H(IJK,'V_MOMENTUM',Xn,Yn,Zn,DELH_n,Nx,Ny,Nz)
660
661 IF((DELH_nc == UNDEFINED).OR.(DELH_n == UNDEFINED)) THEN
662 ALPHA_Vn_c(IJK) = ONE
663 ELSE
664 ALPHA_Vn_c(IJK) = DMIN1(ALPHA_MAX , DELH_nc / DELH_n )
665 ENDIF
666
667 IF(BLOCKED_V_CELL_AT(IJPK)) ALPHA_Vn_c(IJK) = ZERO
668 IF(WALL_V_AT(IJK).AND.WALL_V_AT(IJPK)) ALPHA_Vn_c(IJK) = ZERO
669 IF(J == JEND1) ALPHA_Vn_c(IJK) = ONE
670
671 IF(ALPHA_Vn_c(IJK)<ZERO) THEN
672 WRITE(*,*) 'NEGATIVE ALPHA_Vn_c at IJK=',IJK
673 WRITE(*,*) 'MFIX WILL EXIT NOW.'
674 CALL MFIX_EXIT(MYPE)
675 ENDIF
676
677
678
679
680
681 = X_V(IJPK) - X_V(IJK)
682 Sy = Y_V(IJPK) - Y_V(IJK)
683 Sz = Z_V(IJPK) - Z_V(IJK)
684
685 NOC_V_N(IJK) = (Sx * Nx + Sz * Nz)/(Sy * DELH_n)
686
687 IF(BLOCKED_V_CELL_AT(IJPK)) NOC_V_N(IJK) = ZERO
688 IF(WALL_V_AT(IJK).AND.WALL_V_AT(IJPK)) NOC_V_N(IJK) = ZERO
689 IF(J == JEND1) NOC_V_N(IJK) = ZERO
690
691 IF(DO_K) THEN
692
693
694
695
696
697 (IJK) = DELZ_Vt(IJK) / (DELZ_Vt(IJK) + DELZ_Vb(IJKP))
698 Theta_Vt_bar(IJK) = ONE - Theta_Vt(IJK)
699
700 IF ((Theta_Vt(IJK)>=ONE).OR.(Theta_Vt(IJK)<=ZERO)) THEN
701 Theta_Vt(IJK) = HALF
702 Theta_Vt_bar(IJK) = HALF
703 ENDIF
704
705
706
707
708 = Theta_Vt_bar(IJK) * X_V(IJK) + Theta_Vt(IJK) * X_V(IJKP)
709 Yt = Theta_Vt_bar(IJK) * Y_V(IJK) + Theta_Vt(IJK) * Y_V(IJKP)
710 Zt = Z_NODE(8)
711
712 CALL GET_DEL_H(IJK,'V_MOMENTUM',X_V_tc(IJK),Y_V_tc(IJK),Z_V_tc(IJK),DELH_tc,Nx,Ny,Nz)
713
714 CALL GET_DEL_H(IJK,'V_MOMENTUM',Xt,Yt,Zt,DELH_t,Nx,Ny,Nz)
715
716 IF((DELH_tc == UNDEFINED).OR.(DELH_t == UNDEFINED)) THEN
717 ALPHA_Vt_c(IJK) = ONE
718 ELSE
719 ALPHA_Vt_c(IJK) = DMIN1(ALPHA_MAX , DELH_tc / DELH_t )
720 ENDIF
721
722 IF(BLOCKED_V_CELL_AT(IJKP)) ALPHA_Vt_c(IJK) = ZERO
723 IF(WALL_V_AT(IJK).AND.WALL_V_AT(IJKP)) ALPHA_Vt_c(IJK) = ZERO
724 IF(K == KEND1) ALPHA_Vt_c(IJK) = ONE
725
726 IF(ALPHA_Vt_c(IJK)<ZERO) ALPHA_Vt_c(IJK) = ONE
727
728
729
730
731
732 = X_V(IJKP) - X_V(IJK)
733 Sy = Y_V(IJKP) - Y_V(IJK)
734 Sz = Z_V(IJKP) - Z_V(IJK)
735
736 NOC_V_T(IJK) = (Sx * Nx + Sy * Ny)/(Sz * DELH_t)
737
738
739 IF(BLOCKED_V_CELL_AT(IJKP)) NOC_V_T(IJK) = ZERO
740 IF(WALL_V_AT(IJK).AND.WALL_V_AT(IJKP)) NOC_V_T(IJK) = ZERO
741 IF(K == KEND1) NOC_V_T(IJK) = ZERO
742
743 ENDIF
744
745 ENDIF
746
747
748
749
750
751 IF ( CUT_V_CELL_AT(IJK) ) THEN
752
753 SELECT CASE (PG_OPTION)
754
755 CASE(0)
756
757
758 CASE(1)
759
760 A_VPG_N(IJK) = dmax1(AXZ_V(IJK),AXZ_V(IJMK))
761 A_VPG_S(IJK) = A_VPG_N(IJK)
762
763 CASE(2)
764
765 A_VPG_N(IJK) = AXZ_V(IJK)
766 A_VPG_S(IJK) = AXZ_V(IJMK)
767
768
769 CASE DEFAULT
770
771 WRITE(*,*)'INVALID PRESSURE GRADIENT OPTION:',PG_OPTION
772 WRITE(*,*)'PG_OPTION SHOULD BE SET EQUAL TO 0, 1 OR 2.'
773 WRITE(*,*)'PLEASE VERIFY MFIX.DAT FILE AND TRY AGAIN.'
774 WRITE(*,*) 'MFIX WILL EXIT NOW.'
775 CALL MFIX_EXIT(myPE)
776
777 END SELECT
778
779
780 IF (BLOCKED_CELL_AT(IJK).OR.BLOCKED_CELL_AT(IJPK)) THEN
781 IF(.NOT.WALL_V_AT(IJK)) THEN
782 WALL_V_AT(IJK) = .TRUE.
783 FLAG_N(IJK) = 0
784 IF(PRINT_WARNINGS) THEN
785 WRITE(*,*) 'WARNING: ONLY ONE PRESSURE NODE DETECTED IN V-MOMENTUM CELL IJK =',IJK
786 WRITE(*,*) 'RESETTING U-MOMENTUM CELL AS V_WALL CELL.'
787 ENDIF
788
789
790 ENDIF
791 ENDIF
792
793 ELSE
794
795 A_VPG_N(IJK) = AXZ_V(IJK)
796 A_VPG_S(IJK) = AXZ_V(IJK)
797
798 ENDIF
799
800 ENDIF
801 END DO
802
803 IF(PG_OPTION==0) THEN
804 A_VPG_N = AXZ
805 A_VPG_S = AXZ
806 ENDIF
807
808 RETURN
809
810 END SUBROUTINE GET_3D_ALPHA_V_CUT_CELL
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826 SUBROUTINE GET_3D_ALPHA_W_CUT_CELL
827
828 USE bc
829 USE compar, ONLY: iend1, jend1, kend1, IJKSTART3, IJKEND3, mype, pe_io
830 USE cutcell
831 USE exit, only: mfix_exit
832 USE functions, ONLY: FUNIJK
833 USE geometry, ONLY: axy, axy_w, flag_t
834 USE indices, ONLY: I_OF, J_OF, K_OF
835 USE param1, ONLY: half, one, undefined, zero
836
837 IMPLICIT NONE
838 DOUBLE PRECISION:: Xe,Ye,Ze,Xn,Yn,Zn,Xt,Yt,Zt
839 INTEGER :: I,J,K,IM,JM,KM,IP,JP,KP,IJK,IMJK,IJMK,IPJK,IJPK,IJKM,IJKP
840 DOUBLE PRECISION :: Zi
841 DOUBLE PRECISION :: Sx,Sy,Sz
842 DOUBLE PRECISION :: Nx,Ny,Nz
843 DOUBLE PRECISION :: DELH_ec , DELH_e
844 DOUBLE PRECISION :: DELH_nc , DELH_n
845 DOUBLE PRECISION :: DELH_tc , DELH_t
846 LOGICAL :: U_NODE_AT_TE, U_NODE_AT_BE
847 LOGICAL :: V_NODE_AT_TN, V_NODE_AT_BN
848 INTEGER :: BCV
849 INTEGER ::BCT
850
851 IF(MyPE == PE_IO) THEN
852 WRITE(*,10)'COMPUTING INTERPOLATION FACTORS IN W-MOMENTUM CELLS...'
853 ENDIF
854 10 FORMAT(1X,A)
855
856 DELH_W = UNDEFINED
857
858 Theta_W_te = HALF
859 Theta_W_be = HALF
860 Theta_W_tn = HALF
861 Theta_W_bn = HALF
862 Theta_Wt = HALF
863 Theta_Wt_bar = HALF
864 Theta_We = HALF
865 Theta_We_bar = HALF
866 ALPHA_We_c = ONE
867 NOC_W_E = ZERO
868 Theta_Wn = HALF
869 Theta_Wn_bar = HALF
870 ALPHA_Wn_c = ONE
871 NOC_W_N = ZERO
872 ALPHA_Wt_c = ONE
873 NOC_W_T = ZERO
874 A_WPG_T = AXY_W
875 A_WPG_B = AXY_W
876
877 DO IJK = IJKSTART3, IJKEND3
878 IF(INTERIOR_CELL_AT(IJK)) THEN
879
880 I = I_OF(IJK)
881 J = J_OF(IJK)
882 K = K_OF(IJK)
883
884 IM = I - 1
885 JM = J - 1
886 KM = K - 1
887
888 IP = I + 1
889 JP = J + 1
890 KP = K + 1
891
892 IMJK = FUNIJK(IM,J,K)
893 IPJK = FUNIJK(IP,J,K)
894 IJMK = FUNIJK(I,JM,K)
895 IJPK = FUNIJK(I,JP,K)
896 IJKM = FUNIJK(I,J,KM)
897 IJKP = FUNIJK(I,J,KP)
898
899 CALL GET_CELL_NODE_COORDINATES(IJK,'W_MOMENTUM')
900
901
902
903
904
905 = ((.NOT.BLOCKED_U_CELL_AT(IJKP)).AND.(.NOT.WALL_U_AT(IJKP)))
906 U_NODE_AT_BE = ((.NOT.BLOCKED_U_CELL_AT(IJK)).AND.(.NOT.WALL_U_AT(IJK)))
907
908
909 IF(U_NODE_AT_TE.AND.U_NODE_AT_BE) THEN
910 Theta_W_te(IJK) = (Z_W_ec(IJK) - Z_U(IJK) ) / (Z_U(IJKP) - Z_U(IJK))
911 Theta_W_be(IJK) = ONE - Theta_W_te(IJK)
912 ELSE IF (U_NODE_AT_BE.AND.(.NOT.U_NODE_AT_TE)) THEN
913 Zi = HALF * (Zt_W_int(IJK) + Zt_W_int(IJMK))
914 Theta_W_te(IJK) = (Z_W_ec(IJK) - Z_U(IJK) ) / (Zi - Z_U(IJK))
915 Theta_W_be(IJK) = ONE - Theta_W_te(IJK)
916 ELSE IF ((.NOT.U_NODE_AT_BE).AND.U_NODE_AT_TE) THEN
917 Zi = HALF * (Zt_W_int(IJK) + Zt_W_int(IJMK))
918 Theta_W_te(IJK) = (Z_W_ec(IJK) - Zi) / (Z_U(IJKP) - Zi)
919 Theta_W_be(IJK) = ONE - Theta_W_te(IJK)
920 ELSE
921 Theta_W_te(IJK) = ZERO
922 Theta_W_be(IJK) = ZERO
923 ENDIF
924
925 IF ((Theta_W_te(IJK)>=ONE).OR.(Theta_W_te(IJK)<=ZERO)) THEN
926 Theta_W_te(IJK) = HALF
927 Theta_W_be(IJK) = HALF
928 ENDIF
929
930
931
932
933
934 = ((.NOT.BLOCKED_V_CELL_AT(IJKP)).AND.(.NOT.WALL_V_AT(IJKP)))
935 V_NODE_AT_BN = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.(.NOT.WALL_V_AT(IJK)))
936
937
938 IF(V_NODE_AT_TN.AND.V_NODE_AT_BN) THEN
939 Theta_W_tn(IJK) = (Z_W_nc(IJK) - Z_V(IJK) ) / (Z_V(IJKP) - Z_V(IJK))
940 Theta_W_bn(IJK) = ONE - Theta_W_tn(IJK)
941 ELSE IF (V_NODE_AT_BN.AND.(.NOT.V_NODE_AT_TN)) THEN
942 Zi = HALF * (Zt_W_int(IJK) + Zt_W_int(IMJK))
943 Theta_W_tn(IJK) = (Z_W_nc(IJK) - Z_V(IJK) ) / (Zi - Z_V(IJK))
944 Theta_W_bn(IJK) = ONE - Theta_W_bn(IJK)
945 ELSE IF ((.NOT.V_NODE_AT_BN).AND.V_NODE_AT_TN) THEN
946 Zi = HALF * (Zt_W_int(IJK) + Zt_W_int(IMJK))
947 Theta_W_tn(IJK) = (Z_W_nc(IJK) - Zi) / (Z_V(IJKP) - Zi)
948 Theta_W_bn(IJK) = ONE - Theta_W_bn(IJK)
949 ELSE
950 Theta_W_tn(IJK) = ZERO
951 Theta_W_bn(IJK) = ZERO
952 ENDIF
953
954 IF ((Theta_W_tn(IJK)>=ONE).OR.(Theta_W_tn(IJK)<=ZERO)) THEN
955 Theta_W_tn(IJK) = HALF
956 Theta_W_bn(IJK) = HALF
957 ENDIF
958
959
960
961
962
963 (IJK) = DELZ_Wt(IJK) / (DELZ_Wt(IJK) + DELZ_Wb(IJKP))
964 Theta_Wt_bar(IJK) = ONE - Theta_Wt(IJK)
965
966 IF ((Theta_Wt(IJK)>=ONE).OR.(Theta_Wt(IJK)<=ZERO)) THEN
967 Theta_Wt(IJK) = HALF
968 Theta_Wt_bar(IJK) = HALF
969 ENDIF
970
971 BCV = BC_W_ID(IJK)
972 IF(BCV>0) THEN
973 BCT = BC_TYPE_ENUM(BCV)
974 ELSE
975 BCT = BLANK
976 ENDIF
977
978 IF(BCT ==CG_NSW.OR.BCT ==CG_PSW) THEN
979
980 CALL GET_DEL_H(IJK,'W_MOMENTUM',X_W(IJK),Y_W(IJK),Z_W(IJK),DELH_W(IJK),Nx,Ny,Nz)
981
982 IF(DELH_W(IJK)<ZERO.AND.(.NOT.WALL_W_AT(IJK))) THEN
983 WRITE(*,*) 'NEGATIVE DELH-W AT XYZ=', X_W(IJK),Y_W(IJK),Z_W(IJK)
984 WRITE(*,*) 'DELH_W=', DELH_W(IJK)
985 WRITE(*,*) 'AXY=', AXY(IJK)
986 WRITE(*,*) 'MFIX WILL EXIT NOW.'
987 CALL MFIX_EXIT(MYPE)
988 ENDIF
989
990
991
992
993
994 (IJK) = DELX_We(IJK) / (DELX_We(IJK) + DELX_Ww(IPJK))
995 Theta_We_bar(IJK) = ONE - Theta_We(IJK)
996
997 IF ((Theta_We(IJK)>=ONE).OR.(Theta_We(IJK)<=ZERO)) THEN
998 Theta_We(IJK) = HALF
999 Theta_We_bar(IJK) = HALF
1000 ENDIF
1001
1002
1003
1004 = X_NODE(8)
1005 Ye = Theta_We_bar(IJK) * Y_W(IJK) +Theta_We(IJK) * Y_W(IPJK)
1006 Ze = Theta_We_bar(IJK) * Z_W(IJK) +Theta_We(IJK) * Z_W(IPJK)
1007
1008 CALL GET_DEL_H(IJK,'W_MOMENTUM',X_W_ec(IJK),Y_W_ec(IJK),Z_W_ec(IJK),DELH_ec,Nx,Ny,Nz)
1009
1010 CALL GET_DEL_H(IJK,'W_MOMENTUM',Xe,Ye,Ze,DELH_e,Nx,Ny,Nz)
1011
1012
1013 IF((DELH_ec == UNDEFINED).OR.(DELH_e == UNDEFINED)) THEN
1014 ALPHA_We_c(IJK) = ONE
1015 ELSE
1016 ALPHA_We_c(IJK) = DMIN1(ALPHA_MAX , DELH_ec / DELH_e )
1017 ENDIF
1018
1019 IF(BLOCKED_W_CELL_AT(IPJK)) ALPHA_We_c(IJK) = ZERO
1020 IF(WALL_W_AT(IJK).AND.WALL_W_AT(IPJK)) ALPHA_We_c(IJK) = ZERO
1021 IF(I == IEND1) ALPHA_We_c(IJK) = ONE
1022
1023 IF(ALPHA_We_c(IJK)<ZERO) ALPHA_We_c(IJK) = ONE
1024
1025
1026
1027
1028
1029 = X_W(IPJK) - X_W(IJK)
1030 Sy = Y_W(IPJK) - Y_W(IJK)
1031 Sz = Z_W(IPJK) - Z_W(IJK)
1032
1033 NOC_W_E(IJK) = (Sy * Ny + Sz * Nz)/(Sx * DELH_e)
1034
1035
1036 IF(BLOCKED_W_CELL_AT(IPJK)) NOC_W_E(IJK) = ZERO
1037 IF(WALL_W_AT(IJK).AND.WALL_W_AT(IPJK)) NOC_W_E(IJK) = ZERO
1038 IF(I == IEND1) NOC_W_E(IJK) = ZERO
1039
1040
1041
1042
1043
1044 (IJK) = DELY_Wn(IJK) / (DELY_Wn(IJK) + DELY_Ws(IJPK))
1045 Theta_Wn_bar(IJK) = ONE - Theta_Wn(IJK)
1046
1047 IF ((Theta_Wn(IJK)>=ONE).OR.(Theta_Wn(IJK)<=ZERO)) THEN
1048 Theta_Wn(IJK) = HALF
1049 Theta_Wn_bar(IJK) = HALF
1050 ENDIF
1051
1052
1053
1054
1055
1056 = Theta_Wn_bar(IJK) * X_W(IJK) + Theta_Wn(IJK) * X_W(IJPK)
1057 Yn = Y_NODE(8)
1058 Zn = Theta_Wn_bar(IJK) * Z_W(IJK) + Theta_Wn(IJK) * Z_W(IJPK)
1059
1060 CALL GET_DEL_H(IJK,'W_MOMENTUM',X_W_nc(IJK),Y_W_nc(IJK),Z_W_nc(IJK),DELH_nc,Nx,Ny,Nz)
1061
1062 CALL GET_DEL_H(IJK,'W_MOMENTUM',Xn,Yn,Zn,DELH_n,Nx,Ny,Nz)
1063
1064 IF((DELH_nc == UNDEFINED).OR.(DELH_n == UNDEFINED)) THEN
1065 ALPHA_Wn_c(IJK) = ONE
1066 ELSE
1067 ALPHA_Wn_c(IJK) = DMIN1(ALPHA_MAX , DELH_nc / DELH_n )
1068 ENDIF
1069
1070 IF(BLOCKED_W_CELL_AT(IJPK)) ALPHA_Wn_c(IJK) = ZERO
1071 IF(WALL_W_AT(IJK).AND.WALL_W_AT(IJPK)) ALPHA_Wn_c(IJK) = ZERO
1072 IF(J == JEND1) ALPHA_Wn_c(IJK) = ONE
1073
1074 IF(ALPHA_Wn_c(IJK)<ZERO) ALPHA_Wn_c(IJK) = ONE
1075
1076
1077
1078
1079
1080 = X_W(IJPK) - X_W(IJK)
1081 Sy = Y_W(IJPK) - Y_W(IJK)
1082 Sz = Z_W(IJPK) - Z_W(IJK)
1083
1084 NOC_W_N(IJK) = (Sx * Nx + Sz * Nz)/(Sy * DELH_n)
1085
1086 IF(BLOCKED_W_CELL_AT(IJPK)) NOC_W_N(IJK) = ZERO
1087 IF(WALL_W_AT(IJK).AND.WALL_W_AT(IJPK)) NOC_W_N(IJK) = ZERO
1088 IF(J == JEND1) NOC_W_N(IJK) = ZERO
1089
1090
1091
1092
1093
1094
1095
1096
1097 = Theta_Wt_bar(IJK) * X_W(IJK) + Theta_Wt(IJK) * X_W(IJKP)
1098 Yt = Theta_Wt_bar(IJK) * Y_W(IJK) + Theta_Wt(IJK) * Y_W(IJKP)
1099 Zt = Z_NODE(8)
1100
1101 CALL GET_DEL_H(IJK,'W_MOMENTUM',X_W_tc(IJK),Y_W_tc(IJK),Z_W_tc(IJK),DELH_tc,Nx,Ny,Nz)
1102
1103 CALL GET_DEL_H(IJK,'W_MOMENTUM',Xt,Yt,Zt,DELH_t,Nx,Ny,Nz)
1104
1105 IF((DELH_tc == UNDEFINED).OR.(DELH_t == UNDEFINED)) THEN
1106 ALPHA_Wt_c(IJK) = ONE
1107 ELSE
1108 ALPHA_Wt_c(IJK) = DMIN1(ALPHA_MAX , DELH_tc / DELH_t )
1109 ENDIF
1110
1111 IF(BLOCKED_W_CELL_AT(IJKP)) ALPHA_Wt_c(IJK) = ZERO
1112 IF(WALL_W_AT(IJK).AND.WALL_W_AT(IJKP)) ALPHA_Wt_c(IJK) = ZERO
1113 IF(K == KEND1) ALPHA_Wt_c(IJK) = ONE
1114
1115 IF(ALPHA_Wt_c(IJK)<ZERO) THEN
1116 WRITE(*,*) 'NEGATIVE ALPHA_Wt_c at IJK=',IJK
1117 WRITE(*,*) 'MFIX WILL EXIT NOW.'
1118 CALL MFIX_EXIT(MYPE)
1119 ENDIF
1120
1121
1122
1123
1124
1125 = X_W(IJKP) - X_W(IJK)
1126 Sy = Y_W(IJKP) - Y_W(IJK)
1127 Sz = Z_W(IJKP) - Z_W(IJK)
1128
1129 NOC_W_T(IJK) = (Sx * Nx + Sy * Ny)/(Sz * DELH_t)
1130
1131 IF(BLOCKED_W_CELL_AT(IJKP)) NOC_W_T(IJK) = ZERO
1132 IF(WALL_W_AT(IJK).AND.WALL_W_AT(IJKP)) NOC_W_T(IJK) = ZERO
1133 IF(K == KEND1) NOC_W_T(IJK) = ZERO
1134
1135 ENDIF
1136
1137
1138
1139
1140
1141 IF ( CUT_W_CELL_AT(IJK) ) THEN
1142
1143 SELECT CASE (PG_OPTION)
1144
1145 CASE(0)
1146
1147
1148 CASE(1)
1149
1150 A_WPG_T(IJK) = dmax1(AXY_W(IJK),AXY_W(IJKM))
1151 A_WPG_B(IJK) = A_WPG_T(IJK)
1152
1153 CASE(2)
1154
1155 A_WPG_T(IJK) = AXY_W(IJK)
1156 A_WPG_B(IJK) = AXY_W(IJKM)
1157
1158 CASE DEFAULT
1159
1160 WRITE(*,*)'INVALID PRESSURE GRADIENT OPTION:',PG_OPTION
1161 WRITE(*,*)'PG_OPTION SHOULD BE SET EQUAL TO 0, 1 OR 2.'
1162 WRITE(*,*)'PLEASE VERIFY MFIX.DAT FILE AND TRY AGAIN.'
1163 WRITE(*,*) 'MFIX WILL EXIT NOW.'
1164 CALL MFIX_EXIT(myPE)
1165
1166 END SELECT
1167
1168 IF (BLOCKED_CELL_AT(IJK).OR.BLOCKED_CELL_AT(IJKP)) THEN
1169 IF(.NOT.WALL_W_AT(IJK)) THEN
1170 WALL_W_AT(IJK) = .TRUE.
1171 FLAG_T(IJK) = 0
1172 IF(PRINT_WARNINGS) THEN
1173 WRITE(*,*) 'WARNING: ONLY ONE PRESSURE NODE DETECTED IN W-MOMENTUM CELL IJK =',IJK
1174 WRITE(*,*) 'RESETTING U-MOMENTUM CELL AS W_WALL CELL.'
1175 ENDIF
1176
1177
1178 ENDIF
1179 ENDIF
1180
1181 ELSE
1182
1183 A_WPG_T(IJK) = AXY_W(IJK)
1184 A_WPG_B(IJK) = AXY_W(IJK)
1185
1186 ENDIF
1187
1188 ENDIF
1189 END DO
1190
1191 IF(PG_OPTION==0) THEN
1192 A_WPG_T = AXY
1193 A_WPG_B = AXY
1194 ENDIF
1195
1196 RETURN
1197
1198 END SUBROUTINE GET_3D_ALPHA_W_CUT_CELL
1199