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