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