File: /nfs/home/0/users/jenkins/mfix.git/model/tau_w_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 SUBROUTINE CALC_TAU_W_G(TAU_W_G)
24
25
26
27
28
29
30
31
32 USE param
33 USE param1
34 USE parallel
35 USE matrix
36 USE scales
37 USE constant
38 USE physprop
39 USE fldvar
40 USE visc_g
41 USE rxns
42 USE run
43 USE toleranc
44 USE geometry
45 USE indices
46 USE is
47 USE sendrecv
48 USE compar
49 USE fun_avg
50 USE functions
51
52
53
54
55 USE bc
56 USE quadric
57 USE cutcell
58
59
60
61 IMPLICIT NONE
62
63
64
65
66
67
68
69
70
71 INTEGER IER
72
73
74 DOUBLE PRECISION TAU_W_g(DIMENSION_3)
75
76
77 INTEGER IJK, J, I, IM, IJKP, IMJK, IJKN, IJKNT, IJKS,&
78 IJKST, IJMKP, IJMK, IJKE, IJKTE, IJKW, IJKTW,&
79 IMJKP, K, IJKT, JM, KP, IJKM
80
81
82 DOUBLE PRECISION EPGA
83
84
85 DOUBLE PRECISION duodz
86
87
88 DOUBLE PRECISION Sbv, Ssx, Ssy, Ssz
89
90
91 DOUBLE PRECISION Vxz
92
93
94
95
96 DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
97 LOGICAL :: U_NODE_AT_ET,U_NODE_AT_EB,U_NODE_AT_WT,U_NODE_AT_WB
98 LOGICAL :: V_NODE_AT_NT,V_NODE_AT_NB,V_NODE_AT_ST,V_NODE_AT_SB
99
100 DOUBLE PRECISION :: dudz_at_E,dudz_at_W
101 DOUBLE PRECISION :: dvdz_at_N,dvdz_at_S
102 DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Vi,Sx,Sy,Sz
103 DOUBLE PRECISION :: MU_GT_CUT,SSX_CUT,SSY_CUT
104 DOUBLE PRECISION :: UW_g,VW_g,WW_g
105 INTEGER :: BCV
106 CHARACTER(LEN=9) :: BCT
107
108
109
110
111
112
113
114
115
116
117
118 DO IJK = IJKSTART3, IJKEND3
119 K = K_OF(IJK)
120 IJKT = TOP_OF(IJK)
121 EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
122 IF ( .NOT.IP_AT_T(IJK) .AND. EPGA>DIL_EP_S) THEN
123 J = J_OF(IJK)
124 I = I_OF(IJK)
125 IM = IM1(I)
126 JM = JM1(J)
127 IJKP = KP_OF(IJK)
128 IMJK = IM_OF(IJK)
129 IJKN = NORTH_OF(IJK)
130 IJKNT = TOP_OF(IJKN)
131 IJKS = SOUTH_OF(IJK)
132 IJKST = TOP_OF(IJKS)
133 IJMKP = JM_OF(IJKP)
134 IJMK = JM_OF(IJK)
135 IJKE = EAST_OF(IJK)
136 IJKTE = EAST_OF(IJKT)
137 IJKW = WEST_OF(IJK)
138 IJKTW = WEST_OF(IJKT)
139 IMJKP = KP_OF(IMJK)
140 KP = KP1(K)
141 IJKM = KM_OF(IJK)
142
143
144
145
146
147 IF((.NOT.CARTESIAN_GRID).OR.(CG_SAFE_MODE(5)==1)) THEN
148
149
150
151
152 = (LAMBDA_GT(IJKT)*TRD_G(IJKT)-LAMBDA_GT(IJK)*TRD_G(IJK))*AXY(&
153 IJK)
154
155
156 = AVG_Z_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT)&
157 ,MU_GT(IJKTE),I),K)*(U_G(IJKP)-U_G(IJK))*OX_E(I)*ODZ_T(K)*AYZ_W(&
158 IJK) - AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(&
159 IJKTW),MU_GT(IJKT),IM),K)*(U_G(IMJKP)-U_G(IMJK))*ODZ_T(K)*DY(J)*&
160 (HALF*(DZ(K)+DZ(KP)))
161
162 = AVG_Z_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT)&
163 ,MU_GT(IJKNT),J),K)*(V_G(IJKP)-V_G(IJK))*OX(I)*ODZ_T(K)*AXZ_W(&
164 IJK) - AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(&
165 IJKST),MU_GT(IJKT),JM),K)*(V_G(IJMKP)-V_G(IJMK))*OX(I)*ODZ_T(K)*&
166 AXZ_W(IJMK)
167
168 = MU_GT(IJKT)*(W_G(IJKP)-W_G(IJK))*OX(I)*ODZ(KP)*AXY_W(IJK) - &
169 MU_GT(IJK)*(W_G(IJK)-W_G(IJKM))*OX(I)*ODZ(K)*AXY_W(IJKM)
170
171 ELSE
172
173
174
175
176 = (LAMBDA_GT(IJKT)*TRD_G(IJKT)) * AXY_W(IJK) &
177 -(LAMBDA_GT(IJK) *TRD_G(IJK) ) * AXY_W(IJKM)
178
179
180
181 IF(.NOT.CUT_W_CELL_AT(IJK)) THEN
182
183 SSX = AVG_Z_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT)&
184 ,MU_GT(IJKTE),I),K)*(U_G(IJKP)-U_G(IJK))*ONEoDZ_T_U(IJK)*AYZ_W(&
185 IJK) - AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(&
186 IJKTW),MU_GT(IJKT),IM),K)*(U_G(IMJKP)-U_G(IMJK))*ONEoDZ_T_U(IMJK)*AYZ_W(IMJK)
187
188 SSY = AVG_Z_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT)&
189 ,MU_GT(IJKNT),J),K)*(V_G(IJKP)-V_G(IJK))*ONEoDZ_T_V(IJK)*AXZ_W(&
190 IJK) - AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(&
191 IJKST),MU_GT(IJKT),JM),K)*(V_G(IJMKP)-V_G(IJMK))*ONEoDZ_T_V(IJMK)*&
192 AXZ_W(IJMK)
193
194 = MU_GT(IJKT)*(W_G(IJKP)-W_G(IJK))*ONEoDZ_T_W(IJK)*AXY_W(IJK) - &
195 MU_GT(IJK)*(W_G(IJK)-W_G(IJKM))*ONEoDZ_T_W(IJKM)*AXY_W(IJKM)
196
197 ELSE
198
199 = BC_W_ID(IJK)
200
201 IF(BCV > 0 ) THEN
202 BCT = BC_TYPE(BCV)
203 ELSE
204 BCT = 'NONE'
205 ENDIF
206
207 SELECT CASE (BCT)
208 CASE ('CG_NSW')
209 CUT_TAU_WG = .TRUE.
210 NOC_WG = .TRUE.
211 UW_g = ZERO
212 VW_g = ZERO
213 WW_g = ZERO
214 CASE ('CG_FSW')
215 CUT_TAU_WG = .FALSE.
216 NOC_WG = .FALSE.
217 UW_g = ZERO
218 VW_g = ZERO
219 WW_g = ZERO
220 CASE('CG_PSW')
221 IF(BC_HW_G(BC_W_ID(IJK))==UNDEFINED) THEN
222 = .TRUE.
223 NOC_WG = .TRUE.
224 UW_g = BC_UW_G(BCV)
225 VW_g = BC_VW_G(BCV)
226 WW_g = BC_WW_G(BCV)
227 ELSEIF(BC_HW_G(BC_W_ID(IJK))==ZERO) THEN
228 = .FALSE.
229 NOC_WG = .FALSE.
230 UW_g = ZERO
231 VW_g = ZERO
232 WW_g = ZERO
233 ELSE
234 = .FALSE.
235 NOC_WG = .FALSE.
236 ENDIF
237 CASE ('NONE')
238 TAU_W_G(IJK) = ZERO
239 CYCLE
240 END SELECT
241
242 IF(CUT_TAU_WG) THEN
243 MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJKP)*MU_GT(IJKT))/(VOL(IJK) + VOL(IJKP))
244 ELSE
245 MU_GT_CUT = ZERO
246 ENDIF
247
248
249 = ((.NOT.BLOCKED_U_CELL_AT(IJKP)).AND.(.NOT.WALL_U_AT(IJKP)))
250 U_NODE_AT_EB = ((.NOT.BLOCKED_U_CELL_AT(IJK)).AND.(.NOT.WALL_U_AT(IJK)))
251 U_NODE_AT_WT = ((.NOT.BLOCKED_U_CELL_AT(IMJKP)).AND.(.NOT.WALL_U_AT(IMJKP)))
252 U_NODE_AT_WB = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.(.NOT.WALL_U_AT(IMJK)))
253
254 IF(U_NODE_AT_ET.AND.U_NODE_AT_EB) THEN
255
256 Ui = HALF * (U_G(IJKP) + U_G(IJK))
257 Xi = HALF * (X_U(IJKP) + X_U(IJK))
258 Yi = HALF * (Y_U(IJKP) + Y_U(IJK))
259 Zi = HALF * (Z_U(IJKP) + Z_U(IJK))
260 Sx = X_U(IJKP) - X_U(IJK)
261 Sy = Y_U(IJKP) - Y_U(IJK)
262 Sz = Z_U(IJKP) - Z_U(IJK)
263
264 CALL GET_DEL_H(IJK,'W_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
265
266 dudz_at_E = (U_G(IJKP) - U_G(IJK)) * ONEoDZ_T_U(IJK)
267
268 IF(NOC_WG) dudz_at_E = dudz_at_E - ((Ui-UW_g) * ONEoDZ_T_U(IJK)/DEL_H*(Sx*Nx+Sy*Ny))
269
270 ELSE
271 dudz_at_E = ZERO
272 ENDIF
273
274 IF(U_NODE_AT_WT.AND.U_NODE_AT_WB) THEN
275
276 Ui = HALF * (U_G(IMJKP) + U_G(IMJK))
277 Xi = HALF * (X_U(IMJKP) + X_U(IMJK))
278 Yi = HALF * (Y_U(IMJKP) + Y_U(IMJK))
279 Zi = HALF * (Z_U(IMJKP) + Z_U(IMJK))
280 Sx = X_U(IMJKP) - X_U(IMJK)
281 Sy = Y_U(IMJKP) - Y_U(IMJK)
282 Sz = Z_U(IMJKP) - Z_U(IMJK)
283
284 CALL GET_DEL_H(IJK,'W_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
285
286 dudz_at_W = (U_G(IMJKP) - U_G(IMJK)) * ONEoDZ_T_U(IMJK)
287
288 IF(NOC_WG) dudz_at_W = dudz_at_W - ((Ui-UW_g) * ONEoDZ_T_U(IMJK)/DEL_H*(Sx*Nx+Sy*Ny))
289
290 ELSE
291 dudz_at_W = ZERO
292 ENDIF
293
294 IF(U_NODE_AT_EB) THEN
295 CALL GET_DEL_H(IJK,'W_MOMENTUM',X_U(IJK),Y_U(IJK),Z_U(IJK),Del_H,Nx,Ny,Nz)
296 SSX_CUT = - MU_GT_CUT * (U_G(IJK) - UW_g) / DEL_H * (Nz*Nx) * Area_W_CUT(IJK)
297 ELSE
298 SSX_CUT = ZERO
299 ENDIF
300
301 SSX = AVG_Z_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT)&
302 ,MU_GT(IJKTE),I),K)*dudz_at_E*AYZ_W(&
303 IJK) - AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(&
304 IJKTW),MU_GT(IJKT),IM),K)*dudz_at_W*AYZ_W(IMJK) &
305 + SSX_CUT
306
307
308
309 = ((.NOT.BLOCKED_V_CELL_AT(IJKP)).AND.(.NOT.WALL_V_AT(IJKP)))
310 V_NODE_AT_NB = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.(.NOT.WALL_V_AT(IJK)))
311 V_NODE_AT_ST = ((.NOT.BLOCKED_V_CELL_AT(IJMKP)).AND.(.NOT.WALL_V_AT(IJMKP)))
312 V_NODE_AT_SB = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.(.NOT.WALL_V_AT(IJMK)))
313
314 IF(V_NODE_AT_NT.AND.V_NODE_AT_NB) THEN
315
316 Vi = HALF * (V_G(IJKP) + V_G(IJK))
317 Xi = HALF * (X_V(IJKP) + X_V(IJK))
318 Yi = HALF * (Y_V(IJKP) + Y_V(IJK))
319 Zi = HALF * (Z_V(IJKP) + Z_V(IJK))
320 Sx = X_V(IJKP) - X_V(IJK)
321 Sy = Y_V(IJKP) - Y_V(IJK)
322 Sz = Z_V(IJKP) - Z_V(IJK)
323
324
325 CALL GET_DEL_H(IJK,'W_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
326
327 dvdz_at_N = (V_G(IJKP) - V_G(IJK)) * ONEoDZ_T_V(IJK)
328
329 IF(NOC_WG) dvdz_at_N = dvdz_at_N - ((Vi-VW_g) * ONEoDZ_T_V(IJK)/DEL_H*(Sx*Nx+Sy*Ny))
330
331 ELSE
332 dvdz_at_N = ZERO
333 ENDIF
334
335 IF(V_NODE_AT_ST.AND.V_NODE_AT_SB) THEN
336
337 Vi = HALF * (V_G(IJMKP) + V_G(IJMK))
338 Xi = HALF * (X_V(IJMKP) + X_V(IJMK))
339 Yi = HALF * (Y_V(IJMKP) + Y_V(IJMK))
340 Zi = HALF * (Z_V(IJMKP) + Z_V(IJMK))
341 Sx = X_V(IJMKP) - X_V(IJMK)
342 Sy = Y_V(IJMKP) - Y_V(IJMK)
343 Sz = Z_V(IJMKP) - Z_V(IJMK)
344
345 CALL GET_DEL_H(IJK,'W_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
346
347 dvdz_at_S = (V_G(IJMKP) - V_G(IJMK)) * ONEoDZ_T_V(IJMK)
348
349 IF(NOC_WG) dvdz_at_S = dvdz_at_S - ((Vi-VW_g) * ONEoDZ_T_V(IJMK)/DEL_H*(Sx*Nx+Sy*Ny))
350
351 ELSE
352 dvdz_at_S = ZERO
353 ENDIF
354
355 IF(V_NODE_AT_NB) THEN
356 CALL GET_DEL_H(IJK,'W_MOMENTUM',X_V(IJK),Y_V(IJK),Z_V(IJK),Del_H,Nx,Ny,Nz)
357 SSY_CUT = - MU_GT_CUT * (V_G(IJK) - VW_g) / DEL_H * (Nz*Ny) * Area_W_CUT(IJK)
358 ELSE
359 SSY_CUT = ZERO
360 ENDIF
361
362 SSY = AVG_Z_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT)&
363 ,MU_GT(IJKNT),J),K)*dvdz_at_N*AXZ_W(&
364 IJK) - AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(&
365 IJKST),MU_GT(IJKT),JM),K)*dvdz_at_S*AXZ_W(IJMK) &
366 + SSY_CUT
367
368
369 CALL GET_DEL_H(IJK,'W_MOMENTUM',X_W(IJK),Y_W(IJK),Z_W(IJK),Del_H,Nx,Ny,Nz)
370
371 SSZ = MU_GT(IJKT)*(W_G(IJKP)-W_G(IJK))*ONEoDZ_T_W(IJK)*AXY_W(IJK) - &
372 MU_GT(IJK)*(W_G(IJK)-W_G(IJKM))*OX(I)*ONEoDZ_T_W(IJKM)*AXY_W(IJKM) &
373 - MU_GT_CUT * (W_g(IJK) - WW_g) / DEL_H * (Nz**2) * Area_W_CUT(IJK)
374
375
376
377 END IF
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396 ENDIF
397
398
399
400
401
402
403 IF (CYLINDRICAL) THEN
404
405
406 = SSZ + MU_GT(IJKT)*(U_G(IJKP)+U_G(IMJKP))*OX(I)*AXY_W(IJK)&
407 - MU_GT(IJK)*(U_G(IJK)+U_G(IMJK))*OX(I)*AXY_W(IJKM)
408
409
410 IF (OX_E(IM) /= UNDEFINED) THEN
411 DUODZ = (U_G(IMJKP)-U_G(IMJK))*OX_E(IM)*ODZ_T(K)
412 ELSE
413 DUODZ = ZERO
414 ENDIF
415 VXZ = AVG_Z(MU_GT(IJK),MU_GT(IJKT),K)*OX(I)*HALF*((U_G(IJKP)-U_G&
416 (IJK))*OX_E(I)*ODZ_T(K)+DUODZ)
417
418 ELSE
419 VXZ = ZERO
420 ENDIF
421
422
423 (IJK) = SBV + SSX + SSY + SSZ + VXZ*VOL_W(IJK)
424 ELSE
425 TAU_W_G(IJK) = ZERO
426 ENDIF
427 END DO
428 call send_recv(tau_w_g,2)
429 RETURN
430 END SUBROUTINE CALC_TAU_W_G
431
432
433
434
435
436