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