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