File: /nfs/home/0/users/jenkins/mfix.git/model/tau_v_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_V_G(TAU_V_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
62 IMPLICIT NONE
63
64
65
66
67
68
69
70
71 DOUBLE PRECISION TAU_V_g(DIMENSION_3)
72
73
74 INTEGER I, J, K, IJK, IJKN, JP, IM, KM, IJPK, IJMK,&
75 IJKE, IJKNE, IJKW, IJKNW, IMJPK, IMJK, IJKT,&
76 IJKTN, IJKB, IJKBN, IJKM, IJPKM
77
78
79 DOUBLE PRECISION EPGA
80
81
82 DOUBLE PRECISION Sbv, Ssx, Ssy, Ssz
83
84
85
86
87 DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
88 LOGICAL :: U_NODE_AT_NE,U_NODE_AT_NW,U_NODE_AT_SE,U_NODE_AT_SW
89 LOGICAL :: W_NODE_AT_TN,W_NODE_AT_TS,W_NODE_AT_BN,W_NODE_AT_BS
90 DOUBLE PRECISION :: U_SUM,X_SUM,Y_SUM,Z_SUM
91 DOUBLE PRECISION :: dudy_at_E,dudy_at_W
92 DOUBLE PRECISION :: dwdy_at_T,dwdy_at_B
93 DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Wi,Sx,Sy,Sz
94 DOUBLE PRECISION :: MU_GT_CUT,SSX_CUT,SSZ_CUT
95 DOUBLE PRECISION :: UW_g,VW_g,WW_g
96 INTEGER :: N_SUM
97 INTEGER :: BCV
98 CHARACTER(LEN=9) :: BCT
99
100
101
102
103
104
105
106
107
108
109
110 DO IJK = IJKSTART3, IJKEND3
111 J = J_OF(IJK)
112 IJKN = NORTH_OF(IJK)
113 EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
114 IF ( .NOT.IP_AT_N(IJK) .AND. EPGA>DIL_EP_S) THEN
115 JP = JP1(J)
116 I = I_OF(IJK)
117 IM = IM1(I)
118 K = K_OF(IJK)
119 KM = KM1(K)
120 IJPK = JP_OF(IJK)
121 IJMK = JM_OF(IJK)
122 IJKE = EAST_OF(IJK)
123 IJKNE = EAST_OF(IJKN)
124 IJKW = WEST_OF(IJK)
125 IJKNW = NORTH_OF(IJKW)
126 IMJPK = IM_OF(IJPK)
127 IMJK = IM_OF(IJK)
128 IJKT = TOP_OF(IJK)
129 IJKTN = NORTH_OF(IJKT)
130 IJKB = BOTTOM_OF(IJK)
131 IJKBN = NORTH_OF(IJKB)
132 IJKM = KM_OF(IJK)
133 IJPKM = JP_OF(IJKM)
134
135
136
137 IF((.NOT.CARTESIAN_GRID).OR.(CG_SAFE_MODE(4)==1)) THEN
138
139
140
141
142 = (LAMBDA_GT(IJKN)*TRD_G(IJKN)-LAMBDA_GT(IJK)*TRD_G(IJK))*AXZ(&
143 IJK)
144
145
146 = AVG_Y_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKN)&
147 ,MU_GT(IJKNE),I),J)*(U_G(IJPK)-U_G(IJK))*ODY_N(J)*AYZ_V(IJK) - &
148 AVG_Y_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(IJKNW),&
149 MU_GT(IJKN),IM),J)*(U_G(IMJPK)-U_G(IMJK))*ODY_N(J)*AYZ_V(IMJK)
150 SSY = MU_GT(IJKN)*(V_G(IJPK)-V_G(IJK))*ODY(JP)*AXZ_V(IJK) - MU_GT(&
151 IJK)*(V_G(IJK)-V_G(IJMK))*ODY(J)*AXZ_V(IJMK)
152 SSZ = AVG_Y_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(IJKN)&
153 ,MU_GT(IJKTN),K),J)*(W_G(IJPK)-W_G(IJK))*ODY_N(J)*AXY_V(IJK) - &
154 AVG_Y_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT(IJKBN),&
155 MU_GT(IJKN),KM),J)*(W_G(IJPKM)-W_G(IJKM))*ODY_N(J)*AXY_V(IJKM)
156
157 ELSE
158
159
160
161
162 = (LAMBDA_GT(IJKN)*TRD_G(IJKN)) * AXZ_V(IJK) &
163 -(LAMBDA_GT(IJK) *TRD_G(IJK) ) * AXZ_V(IJMK)
164
165
166
167 IF(.NOT.CUT_V_CELL_AT(IJK)) THEN
168
169 SSX = AVG_Y_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKN)&
170 ,MU_GT(IJKNE),I),J)*(U_G(IJPK)-U_G(IJK))*ONEoDY_N_U(IJK)*AYZ_V(IJK) - &
171 AVG_Y_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(IJKNW),&
172 MU_GT(IJKN),IM),J)*(U_G(IMJPK)-U_G(IMJK))*ONEoDY_N_U(IMJK)*AYZ_V(IMJK)
173
174 SSY = MU_GT(IJKN)*(V_G(IJPK)-V_G(IJK))*ONEoDY_N_V(IJK)*AXZ_V(IJK) - MU_GT(&
175 IJK)*(V_G(IJK)-V_G(IJMK))*ONEoDY_N_V(IJMK)*AXZ_V(IJMK)
176
177
178
179 IF(DO_K) THEN
180
181 SSZ = AVG_Y_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(IJKN)&
182 ,MU_GT(IJKTN),K),J)*(W_G(IJPK)-W_G(IJK))*ONEoDY_N_W(IJK)*AXY_V(IJK) - &
183 AVG_Y_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT(IJKBN),&
184 MU_GT(IJKN),KM),J)*(W_G(IJPKM)-W_G(IJKM))*ONEoDY_N_W(IJKM)*AXY_V(IJKM)
185
186 ELSE
187 SSZ = ZERO
188 ENDIF
189
190 ELSE
191
192 = BC_V_ID(IJK)
193
194 IF(BCV > 0 ) THEN
195 BCT = BC_TYPE(BCV)
196 ELSE
197 BCT = 'NONE'
198 ENDIF
199
200 SELECT CASE (BCT)
201 CASE ('CG_NSW')
202 CUT_TAU_VG = .TRUE.
203 NOC_VG = .TRUE.
204 UW_g = ZERO
205 VW_g = ZERO
206 WW_g = ZERO
207 CASE ('CG_FSW')
208 CUT_TAU_VG = .FALSE.
209 NOC_VG = .FALSE.
210 UW_g = ZERO
211 VW_g = ZERO
212 WW_g = ZERO
213 CASE('CG_PSW')
214 IF(BC_HW_G(BC_V_ID(IJK))==UNDEFINED) THEN
215 = .TRUE.
216 NOC_VG = .TRUE.
217 UW_g = BC_UW_G(BCV)
218 VW_g = BC_VW_G(BCV)
219 WW_g = BC_WW_G(BCV)
220 ELSEIF(BC_HW_G(BC_V_ID(IJK))==ZERO) THEN
221 = .FALSE.
222 NOC_VG = .FALSE.
223 UW_g = ZERO
224 VW_g = ZERO
225 WW_g = ZERO
226 ELSE
227 = .FALSE.
228 NOC_VG = .FALSE.
229 ENDIF
230 CASE ('NONE')
231 TAU_V_G(IJK) = ZERO
232 CYCLE
233 END SELECT
234
235
236 IF(CUT_TAU_VG) THEN
237 MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJPK)*MU_GT(IJKN))/(VOL(IJK) + VOL(IJPK))
238 ELSE
239 MU_GT_CUT = ZERO
240 ENDIF
241
242
243
244 = ((.NOT.BLOCKED_U_CELL_AT(IJPK)).AND.(.NOT.WALL_U_AT(IJPK)))
245 U_NODE_AT_SE = ((.NOT.BLOCKED_U_CELL_AT(IJK)).AND.(.NOT.WALL_U_AT(IJK)))
246 U_NODE_AT_NW = ((.NOT.BLOCKED_U_CELL_AT(IMJPK)).AND.(.NOT.WALL_U_AT(IMJPK)))
247 U_NODE_AT_SW = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.(.NOT.WALL_U_AT(IMJK)))
248
249
250 IF(U_NODE_AT_NE.AND.U_NODE_AT_SE) THEN
251
252 Ui = HALF * (U_G(IJPK) + U_G(IJK))
253 Xi = HALF * (X_U(IJPK) + X_U(IJK))
254 Yi = HALF * (Y_U(IJPK) + Y_U(IJK))
255 Zi = HALF * (Z_U(IJPK) + Z_U(IJK))
256 Sx = X_U(IJPK) - X_U(IJK)
257 Sy = Y_U(IJPK) - Y_U(IJK)
258 Sz = Z_U(IJPK) - Z_U(IJK)
259
260
261 CALL GET_DEL_H(IJK,'V_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
262
263 dudy_at_E = (U_G(IJPK) - U_G(IJK)) * ONEoDY_N_U(IJK)
264
265 IF(NOC_VG) dudy_at_E = dudy_at_E - ((Ui-UW_g) * ONEoDY_N_U(IJK)/DEL_H*(Sx*Nx+Sz*Nz))
266
267 ELSE
268 dudy_at_E = ZERO
269 ENDIF
270
271 IF(U_NODE_AT_NW.AND.U_NODE_AT_SW) THEN
272
273 Ui = HALF * (U_G(IMJPK) + U_G(IMJK))
274 Xi = HALF * (X_U(IMJPK) + X_U(IMJK))
275 Yi = HALF * (Y_U(IMJPK) + Y_U(IMJK))
276 Zi = HALF * (Z_U(IMJPK) + Z_U(IMJK))
277 Sx = X_U(IMJPK) - X_U(IMJK)
278 Sy = Y_U(IMJPK) - Y_U(IMJK)
279 Sz = Z_U(IMJPK) - Z_U(IMJK)
280
281 CALL GET_DEL_H(IJK,'V_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
282
283 dudy_at_W = (U_G(IMJPK) - U_G(IMJK)) * ONEoDY_N_U(IMJK)
284
285 IF(NOC_VG) dudy_at_W = dudy_at_W - ((Ui-UW_g) * ONEoDY_N_U(IMJK)/DEL_H*(Sx*Nx+Sz*Nz))
286
287 ELSE
288 dudy_at_W = ZERO
289 ENDIF
290
291 U_SUM = ZERO
292 X_SUM = ZERO
293 Y_SUM = ZERO
294 Z_SUM = ZERO
295 N_SUM = 0
296
297 IF(U_NODE_AT_SE) THEN
298 CALL GET_DEL_H(IJK,'V_MOMENTUM',X_U(IJK),Y_U(IJK),Z_U(IJK),Del_H,Nx,Ny,Nz)
299 SSX_CUT = - MU_GT_CUT * (U_G(IJK) - UW_g) / DEL_H * (Ny*Nx) * Area_V_CUT(IJK)
300 ELSE
301 SSX_CUT = ZERO
302 ENDIF
303
304 SSX = AVG_Y_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKN)&
305 ,MU_GT(IJKNE),I),J)*dudy_at_E*AYZ_V(IJK) - &
306 AVG_Y_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(IJKNW),&
307 MU_GT(IJKN),IM),J)*dudy_at_W*AYZ_V(IMJK) &
308 + SSX_CUT
309
310
311 CALL GET_DEL_H(IJK,'V_MOMENTUM',X_V(IJK),Y_V(IJK),Z_V(IJK),Del_H,Nx,Ny,Nz)
312
313 SSY = MU_GT(IJKN)*(V_G(IJPK)-V_G(IJK))*ONEoDY_N_V(IJK)*AXZ_V(IJK) - MU_GT(&
314 IJK)*(V_G(IJK)-V_G(IJMK))*ONEoDY_N_V(IJMK)*AXZ_V(IJMK) &
315 - MU_GT_CUT * (V_g(IJK) - VW_g) / DEL_H * (Ny**2) * Area_V_CUT(IJK)
316
317
318 IF(DO_K) THEN
319
320 W_NODE_AT_TN = ((.NOT.BLOCKED_W_CELL_AT(IJPK)).AND.(.NOT.WALL_W_AT(IJPK)))
321 W_NODE_AT_TS = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.(.NOT.WALL_W_AT(IJK)))
322 W_NODE_AT_BN = ((.NOT.BLOCKED_W_CELL_AT(IJPKM)).AND.(.NOT.WALL_W_AT(IJPKM)))
323 W_NODE_AT_BS = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.(.NOT.WALL_W_AT(IJKM)))
324
325 IF(W_NODE_AT_TN.AND.W_NODE_AT_TS) THEN
326
327 Wi = HALF * (W_G(IJPK) + W_G(IJK))
328 Xi = HALF * (X_W(IJPK) + X_W(IJK))
329 Yi = HALF * (Y_W(IJPK) + Y_W(IJK))
330 Zi = HALF * (Z_W(IJPK) + Z_W(IJK))
331 Sx = X_W(IJPK) - X_W(IJK)
332 Sy = Y_W(IJPK) - Y_W(IJK)
333 Sz = Z_W(IJPK) - Z_W(IJK)
334
335 CALL GET_DEL_H(IJK,'V_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
336
337 dwdy_at_T = (W_G(IJPK) - W_G(IJK)) * ONEoDY_N_W(IJK)
338
339 IF(NOC_VG) dwdy_at_T = dwdy_at_T - ((Wi-WW_g) * ONEoDY_N_W(IJK)/DEL_H*(Sx*Nx+Sz*Nz))
340
341 ELSE
342 dwdy_at_T = ZERO
343 ENDIF
344
345 IF(W_NODE_AT_BN.AND.W_NODE_AT_BS) THEN
346
347 Wi = HALF * (W_G(IJPKM) + W_G(IJKM))
348 Xi = HALF * (X_W(IJPKM) + X_W(IJKM))
349 Yi = HALF * (Y_W(IJPKM) + Y_W(IJKM))
350 Zi = HALF * (Z_W(IJPKM) + Z_W(IJKM))
351 Sx = X_W(IJPKM) - X_W(IJKM)
352 Sy = Y_W(IJPKM) - Y_W(IJKM)
353 Sz = Z_W(IJPKM) - Z_W(IJKM)
354
355 CALL GET_DEL_H(IJK,'V_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
356
357 dwdy_at_B = (W_G(IJPKM) - W_G(IJKM)) * ONEoDY_N_W(IJKM)
358
359 IF(NOC_VG) dwdy_at_B = dwdy_at_B - ((Wi-WW_g) * ONEoDY_N_W(IJKM)/DEL_H*(Sx*Nx+Sz*Nz))
360
361 ELSE
362 dwdy_at_B = ZERO
363 ENDIF
364
365 IF(W_NODE_AT_TS) THEN
366 CALL GET_DEL_H(IJK,'V_MOMENTUM',X_W(IJK),Y_W(IJK),Z_W(IJK),Del_H,Nx,Ny,Nz)
367 SSZ_CUT = - MU_GT_CUT * (W_G(IJK) - WW_g) / DEL_H * (Ny*Nz) * Area_V_CUT(IJK)
368 ELSE
369 SSZ_CUT = ZERO
370 ENDIF
371
372 SSZ = AVG_Y_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(IJKN)&
373 ,MU_GT(IJKTN),K),J)*dwdy_at_T*AXY_V(IJK) - &
374 AVG_Y_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT(IJKBN),&
375 MU_GT(IJKN),KM),J)*dwdy_at_B*AXY_V(IJKM) &
376 + SSZ_CUT
377
378 ELSE
379
380 SSZ = ZERO
381
382 ENDIF
383
384
385 END IF
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402 ENDIF
403
404
405
406
407
408
409 (IJK) = SBV + SSX + SSY + SSZ
410 ELSE
411 TAU_V_G(IJK) = ZERO
412 ENDIF
413 END DO
414 call send_recv(tau_v_g,2)
415 RETURN
416 END SUBROUTINE CALC_TAU_V_G
417
418
419
420
421
422