File: /nfs/home/0/users/jenkins/mfix.git/model/tau_u_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_U_G(TAU_U_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 compar
48 USE sendrecv
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 DOUBLE PRECISION TAU_U_g(DIMENSION_3)
71
72
73 INTEGER I, J, JM, K, KM, IJK, IJKE, IPJK, IP, IMJK, IJKN,&
74 IJKNE, IJKS, IJKSE, IPJMK, IJMK, IJKT, IJKTE,&
75 IJKB, IJKBE, IJKM, IPJKM
76
77
78 DOUBLE PRECISION EPGA
79
80
81 DOUBLE PRECISION EPMU_gte, EPMU_gbe, EPMUGA
82
83
84 DOUBLE PRECISION dWoXdz
85
86
87 DOUBLE PRECISION Sbv, Ssx, Ssy, Ssz
88
89
90 DOUBLE PRECISION Vtzb
91
92
93
94
95 DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
96 LOGICAL :: V_NODE_AT_NE,V_NODE_AT_NW,V_NODE_AT_SE,V_NODE_AT_SW
97 LOGICAL :: W_NODE_AT_TE,W_NODE_AT_TW,W_NODE_AT_BE,W_NODE_AT_BW
98 DOUBLE PRECISION :: dvdx_at_N,dvdx_at_S
99 DOUBLE PRECISION :: dwdx_at_T,dwdx_at_B
100 DOUBLE PRECISION :: Xi,Yi,Zi,Vi,Wi,Sx,Sy,Sz
101 DOUBLE PRECISION :: MU_GT_CUT,SSY_CUT,SSZ_CUT
102 DOUBLE PRECISION :: UW_g,VW_g,WW_g
103 INTEGER :: BCV
104 CHARACTER(LEN=9) :: BCT
105
106
107
108
109
110
111
112
113
114
115
116 DO IJK = IJKSTART3, IJKEND3
117 I = I_OF(IJK)
118 IJKE = EAST_OF(IJK)
119 EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
120 IF ( .NOT.IP_AT_E(IJK) .AND. EPGA>DIL_EP_S) THEN
121 IP = IP1(I)
122 J = J_OF(IJK)
123 JM = JM1(J)
124 K = K_OF(IJK)
125 KM = KM1(K)
126 IPJK = IP_OF(IJK)
127 IMJK = IM_OF(IJK)
128 IJKN = NORTH_OF(IJK)
129 IJKNE = EAST_OF(IJKN)
130 IJKS = SOUTH_OF(IJK)
131 IJKSE = EAST_OF(IJKS)
132 IPJMK = JM_OF(IPJK)
133 IJMK = JM_OF(IJK)
134 IJKT = TOP_OF(IJK)
135 IJKTE = EAST_OF(IJKT)
136 IJKB = BOTTOM_OF(IJK)
137 IJKBE = EAST_OF(IJKB)
138 IJKM = KM_OF(IJK)
139 IPJKM = IP_OF(IJKM)
140
141
142
143 IF((.NOT.CARTESIAN_GRID).OR.(CG_SAFE_MODE(3)==1)) THEN
144
145
146
147
148 = (LAMBDA_GT(IJKE)*TRD_G(IJKE)-LAMBDA_GT(IJK)*TRD_G(IJK))*AYZ(&
149 IJK)
150
151
152 = MU_GT(IJKE)*(U_G(IPJK)-U_G(IJK))*ODX(IP)*AYZ_U(IJK) - MU_GT(&
153 IJK)*(U_G(IJK)-U_G(IMJK))*ODX(I)*AYZ_U(IMJK)
154 SSY = AVG_X_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKE)&
155 ,MU_GT(IJKNE),J),I)*(V_G(IPJK)-V_G(IJK))*ODX_E(I)*AXZ_U(IJK) - &
156 AVG_X_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(IJKSE),&
157 MU_GT(IJKE),JM),I)*(V_G(IPJMK)-V_G(IJMK))*ODX_E(I)*AXZ_U(IJMK)
158 EPMU_GTE = AVG_X_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(&
159 IJKE),MU_GT(IJKTE),K),I)
160 EPMU_GBE = AVG_X_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT&
161 (IJKBE),MU_GT(IJKE),KM),I)
162 SSZ = EPMU_GTE*(W_G(IPJK)-W_G(IJK))*ODX_E(I)*AXY_U(IJK) - EPMU_GBE*&
163 (W_G(IPJKM)-W_G(IJKM))*ODX_E(I)*AXY_U(IJKM)
164
165 ELSE
166
167
168 = (LAMBDA_GT(IJKE)*TRD_G(IJKE)) * AYZ_U(IJK) &
169 -(LAMBDA_GT(IJK) *TRD_G(IJK) ) * AYZ_U(IMJK)
170
171
172 IF(.NOT.CUT_U_CELL_AT(IJK)) THEN
173
174 SSX = MU_GT(IJKE)*(U_G(IPJK)-U_G(IJK))*ONEoDX_E_U(IJK)*AYZ_U(IJK) - MU_GT(&
175 IJK)*(U_G(IJK)-U_G(IMJK))*ONEoDX_E_U(IMJK)*AYZ_U(IMJK)
176
177 SSY = AVG_X_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKE)&
178 ,MU_GT(IJKNE),J),I)*(V_G(IPJK)-V_G(IJK))*ONEoDX_E_V(IJK)*AXZ_U(IJK) - &
179 AVG_X_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(IJKSE),&
180 MU_GT(IJKE),JM),I)*(V_G(IPJMK)-V_G(IJMK))*ONEoDX_E_V(IJMK)*AXZ_U(IJMK)
181
182 IF(DO_K) THEN
183 EPMU_GTE = AVG_X_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(&
184 IJKE),MU_GT(IJKTE),K),I)
185 EPMU_GBE = AVG_X_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT&
186 (IJKBE),MU_GT(IJKE),KM),I)
187
188 SSZ = EPMU_GTE*(W_G(IPJK)-W_G(IJK))*ONEoDX_E_W(IJK)*AXY_U(IJK) - EPMU_GBE*&
189 (W_G(IPJKM)-W_G(IJKM))*ONEoDX_E_W(IJKM)*AXY_U(IJKM)
190
191 ELSE
192 SSZ = ZERO
193 ENDIF
194
195 ELSE
196
197 BCV = BC_U_ID(IJK)
198
199 IF(BCV > 0 ) THEN
200 BCT = BC_TYPE(BCV)
201 ELSE
202 BCT = 'NONE'
203 ENDIF
204
205 SELECT CASE (BCT)
206
207 CASE ('CG_NSW')
208 CUT_TAU_UG = .TRUE.
209 NOC_UG = .TRUE.
210 UW_g = ZERO
211 VW_g = ZERO
212 WW_g = ZERO
213 CASE ('CG_FSW')
214 CUT_TAU_UG = .FALSE.
215 NOC_UG = .FALSE.
216 UW_g = ZERO
217 VW_g = ZERO
218 WW_g = ZERO
219 CASE('CG_PSW')
220 IF(BC_HW_G(BC_U_ID(IJK))==UNDEFINED) THEN
221 = .TRUE.
222 NOC_UG = .TRUE.
223 UW_g = BC_UW_G(BCV)
224 VW_g = BC_VW_G(BCV)
225 WW_g = BC_WW_G(BCV)
226 ELSEIF(BC_HW_G(BC_U_ID(IJK))==ZERO) THEN
227 = .FALSE.
228 NOC_UG = .FALSE.
229 ELSE
230 = .FALSE.
231 NOC_UG = .FALSE.
232 UW_g = ZERO
233 VW_g = ZERO
234 WW_g = ZERO
235 ENDIF
236 CASE ('NONE')
237 TAU_U_G(IJK) = ZERO
238 CYCLE
239 END SELECT
240
241 IF(CUT_TAU_UG) THEN
242 MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IPJK)*MU_GT(IJKE))/(VOL(IJK) + VOL(IPJK))
243 ELSE
244 MU_GT_CUT = ZERO
245 ENDIF
246
247
248
249 CALL GET_DEL_H(IJK,'U_MOMENTUM',X_U(IJK),Y_U(IJK),Z_U(IJK),Del_H,Nx,Ny,Nz)
250
251 SSX = MU_GT(IJKE)*(U_G(IPJK)-U_G(IJK))*ONEoDX_E_U(IJK)*AYZ_U(IJK) - MU_GT(&
252 IJK)*(U_G(IJK)-U_G(IMJK))*ONEoDX_E_U(IMJK)*AYZ_U(IMJK) &
253 - MU_GT_CUT * (U_g(IJK) - UW_g) / DEL_H * (Nx**2) * Area_U_CUT(IJK)
254
255
256
257 = ((.NOT.BLOCKED_V_CELL_AT(IPJK)).AND.(.NOT.WALL_V_AT(IPJK)))
258 V_NODE_AT_NW = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.(.NOT.WALL_V_AT(IJK)))
259 V_NODE_AT_SE = ((.NOT.BLOCKED_V_CELL_AT(IPJMK)).AND.(.NOT.WALL_V_AT(IPJMK)))
260 V_NODE_AT_SW = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.(.NOT.WALL_V_AT(IJMK)))
261
262 IF(V_NODE_AT_NE.AND.V_NODE_AT_NW) THEN
263
264 Vi = HALF * (V_G(IPJK) + V_G(IJK))
265 Xi = HALF * (X_V(IPJK) + X_V(IJK))
266 Yi = HALF * (Y_V(IPJK) + Y_V(IJK))
267 Zi = HALF * (Z_V(IPJK) + Z_V(IJK))
268 Sx = X_V(IPJK) - X_V(IJK)
269 Sy = Y_V(IPJK) - Y_V(IJK)
270 Sz = Z_V(IPJK) - Z_V(IJK)
271
272 CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
273
274 dvdx_at_N = (V_G(IPJK) - V_G(IJK)) * ONEoDX_E_V(IJK)
275
276 IF(NOC_UG) dvdx_at_N = dvdx_at_N - ((Vi-VW_g)*ONEoDX_E_V(IJK)/DEL_H*(Sy*Ny+Sz*Nz))
277
278 ELSE
279 dvdx_at_N = ZERO
280 ENDIF
281
282 IF(V_NODE_AT_SE.AND.V_NODE_AT_SW) THEN
283
284 Vi = HALF * (V_G(IPJMK) + V_G(IJMK))
285 Xi = HALF * (X_V(IPJMK) + X_V(IJMK))
286 Yi = HALF * (Y_V(IPJMK) + Y_V(IJMK))
287 Zi = HALF * (Z_V(IPJMK) + Z_V(IJMK))
288 Sx = X_V(IPJMK) - X_V(IJMK)
289 Sy = Y_V(IPJMK) - Y_V(IJMK)
290 Sz = Z_V(IPJMK) - Z_V(IJMK)
291
292 CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
293
294 dvdx_at_S = (V_G(IPJMK) - V_G(IJMK)) * ONEoDX_E_V(IJMK)
295
296 IF(NOC_UG) dvdx_at_S = dvdx_at_S - ((Vi-VW_g)*ONEoDX_E_V(IJMK)/DEL_H*(Sy*Ny+Sz*Nz))
297
298 ELSE
299 dvdx_at_S = ZERO
300 ENDIF
301
302 IF(V_NODE_AT_NW) THEN
303 CALL GET_DEL_H(IJK,'U_MOMENTUM',X_V(IJK),Y_V(IJK),Z_V(IJK),Del_H,Nx,Ny,Nz)
304 SSY_CUT = - MU_GT_CUT * (V_G(IJK) - VW_g) / DEL_H * (Nx*Ny) * Area_U_CUT(IJK)
305 ELSE
306 SSY_CUT = ZERO
307 ENDIF
308
309
310 SSY = AVG_X_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKE)&
311 ,MU_GT(IJKNE),J),I)*dvdx_at_N*AXZ_U(IJK) - &
312 AVG_X_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(IJKSE),&
313 MU_GT(IJKE),JM),I)*dvdx_at_S*AXZ_U(IJMK) &
314 + SSY_CUT
315
316
317
318 IF(DO_K) THEN
319
320 W_NODE_AT_TE = ((.NOT.BLOCKED_W_CELL_AT(IPJK)).AND.(.NOT.WALL_W_AT(IPJK)))
321 W_NODE_AT_TW = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.(.NOT.WALL_W_AT(IJK)))
322 W_NODE_AT_BE = ((.NOT.BLOCKED_W_CELL_AT(IPJKM)).AND.(.NOT.WALL_W_AT(IPJKM)))
323 W_NODE_AT_BW = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.(.NOT.WALL_W_AT(IJKM)))
324
325 IF(W_NODE_AT_TE.AND.W_NODE_AT_TW) THEN
326
327 Wi = HALF * (W_G(IPJK) + W_G(IJK))
328 Xi = HALF * (X_W(IPJK) + X_W(IJK))
329 Yi = HALF * (Y_W(IPJK) + Y_W(IJK))
330 Zi = HALF * (Z_W(IPJK) + Z_W(IJK))
331 Sx = X_W(IPJK) - X_W(IJK)
332 Sy = Y_W(IPJK) - Y_W(IJK)
333 Sz = Z_W(IPJK) - Z_W(IJK)
334
335 CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
336
337 dwdx_at_T = (W_G(IPJK) - W_G(IJK)) * ONEoDX_E_W(IJK)
338
339 IF(NOC_UG) dwdx_at_T = dwdx_at_T - ((Wi-WW_g) * ONEoDX_E_W(IJK)/DEL_H*(Sy*Ny+Sz*Nz))
340
341 ELSE
342 dwdx_at_T = ZERO
343 ENDIF
344
345
346 IF(W_NODE_AT_BE.AND.W_NODE_AT_BW) THEN
347
348 Wi = HALF * (W_G(IPJKM) + W_G(IJKM))
349 Xi = HALF * (X_W(IPJKM) + X_W(IJKM))
350 Yi = HALF * (Y_W(IPJKM) + Y_W(IJKM))
351 Zi = HALF * (Z_W(IPJKM) + Z_W(IJKM))
352 Sx = X_W(IPJKM) - X_W(IJKM)
353 Sy = Y_W(IPJKM) - Y_W(IJKM)
354 Sz = Z_W(IPJKM) - Z_W(IJKM)
355
356 CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
357
358 dwdx_at_B = (W_G(IPJKM) - W_G(IJKM)) * ONEoDX_E_W(IJKM)
359
360 IF(NOC_UG) dwdx_at_B = dwdx_at_B - ((Wi-WW_g) * ONEoDX_E_W(IJKM)/DEL_H*(Sy*Ny+Sz*Nz))
361
362 ELSE
363 dwdx_at_B = ZERO
364 ENDIF
365
366 IF(W_NODE_AT_TW) THEN
367 CALL GET_DEL_H(IJK,'U_MOMENTUM',X_W(IJK),Y_W(IJK),Z_W(IJK),Del_H,Nx,Ny,Nz)
368 SSZ_CUT = - MU_GT_CUT * (W_G(IJK) - WW_g) / DEL_H * (Nx*Nz) * Area_U_CUT(IJK)
369 ELSE
370 SSZ_CUT = ZERO
371 ENDIF
372
373 EPMU_GTE = AVG_X_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(&
374 IJKE),MU_GT(IJKTE),K),I)
375 EPMU_GBE = AVG_X_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT&
376 (IJKBE),MU_GT(IJKE),KM),I)
377 SSZ = EPMU_GTE*dwdx_at_T*AXY_U(IJK) &
378 - EPMU_GBE*dwdx_at_B*AXY_U(IJKM) &
379 + SSZ_CUT
380
381 ELSE
382
383 SSZ = ZERO
384
385 ENDIF
386
387 ENDIF
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406 ENDIF
407
408
409
410
411
412 IF (CYLINDRICAL) THEN
413
414
415 = SSZ - (EPMU_GTE*(HALF*(W_G(IPJK)+W_G(IJK))*OX_E(I))*AXY_U(&
416 IJK)-EPMU_GBE*(HALF*(W_G(IPJKM)+W_G(IJKM))*OX_E(I))*AXY_U(&
417 IJKM))
418
419
420 = AVG_X(MU_G(IJK),MU_G(IJKE),I)
421 DWOXDZ = HALF*((W_G(IJK)-W_G(IJKM))*OX(I)*ODZ(K)+(W_G(IPJK)-W_G(&
422 IPJKM))*OX(IP)*ODZ(K))
423 VTZB = -2.d0*EPMUGA*OX_E(I)*DWOXDZ
424 ELSE
425 VTZB = ZERO
426 ENDIF
427
428
429 (IJK) = SBV + SSX + SSY + SSZ + VTZB*VOL_U(IJK)
430 ELSE
431 TAU_U_G(IJK) = ZERO
432 ENDIF
433 END DO
434
435 call send_recv(tau_u_g,2)
436
437 RETURN
438
439 END SUBROUTINE CALC_TAU_U_G
440
441
442
443
444