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