File: /nfs/home/0/users/jenkins/mfix.git/model/tau_v_s.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_Tau_V_s                                            C
4     !  Purpose: Cross terms in the gradient of stress in V_s momentum      c
5     !                                                                      C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 19-DEC-96  C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Revision Number: 1                                                  C
11     !  Purpose: To incorporate Cartesian grid modifications                C
12     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
13     !                                                                      C
14     
15     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
16     
17           SUBROUTINE CALC_TAU_V_S(TAU_V_S)
18     
19     !-----------------------------------------------
20     ! Modules
21     !-----------------------------------------------
22           USE param1, only: zero, half, one
23           USE fun_avg
24     
25     ! number of solids phases
26           USE physprop, only: smax, mmax
27     
28     ! x,y,z-components of solids velocity
29           USE fldvar, only: u_s, v_s, w_s
30     ! solids phase particle bulk density and material density
31           USE fldvar, only: rop_s, ro_s, ep_s
32     
33     ! viscous solids transport coefficients
34           USE visc_s, only: mu_s, lambda_s
35     ! trace of rate of strain tensor
36           USE visc_s, only: trd_s
37     
38     ! dilute threshold
39           USE toleranc, only: dil_ep_s
40     
41     ! kinetic theories
42           USE run, only: kt_type_enum, ghd_2007
43     ! runtime flag for treating the system as if shearing
44           USE run, only: shear, V_sh
45     
46     ! primarily needed for function.inc
47           USE compar
48           USE geometry
49           USE indices
50     
51     ! for sendrecv calls
52           USE sendrecv
53     
54     ! for cutcell:
55     ! wall velocity for slip conditions
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     ! Dummy arguments
66     !-----------------------------------------------
67     ! TAU_V_s
68           DOUBLE PRECISION, INTENT(OUT) :: TAU_V_s(DIMENSION_3, DIMENSION_M)
69     !-----------------------------------------------
70     ! Local Variables
71     !-----------------------------------------------
72     ! Indices
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     ! Phase index
77           INTEGER :: M, L
78     
79     ! Average volume fraction
80           DOUBLE PRECISION :: EPSA, EPStmp
81     
82     ! Source terms (Surface)
83           DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
84     ! Shearing variables
85           DOUBLE PRECISION :: Source_diff,Diffco_e,Diffco_w
86     
87     ! for cartesian grid implementation:
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     !!$omp  parallel do private( IJK, I, IJKE, EPSA, EPStmp,  J,  K, KM,  &
103     !!$omp& JP,IM,IJPK,IJKW,IJKNW,IMJPK,IJKTN,IJKBN,IJPKM, &
104     !!$omp&  IMJK,IJKN,IJKNE,IJMK,IJKT,  &
105     !!$omp&  IJKB,IJKM, &
106     !!$omp&  SBV,  SSX,SSY,   SSZ,&
107     !!$omp&  Source_diff, Diffco_e,Diffco_w)
108     !!$omp&  schedule(static)
109     
110             DO IJK = IJKSTART3, IJKEND3
111     
112     ! Skip walls where some values are undefined.
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     ! NON CARTESIAN GRID CASE
151     ! ---------------------------------------------------------------->>>
152     ! Surface forces
153     
154     ! bulk viscosity term
155                   SBV = (LAMBDA_S(IJKN,M)*TRD_S(IJKN,M)-&
156                          LAMBDA_S(IJK,M)*TRD_S(IJK,M))*AXZ(IJK)
157     
158     ! shear stress terms
159                   SSX = 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     !?              IF(DO_K) THEN
168                   SSZ = 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     !?              ELSE
175     !?                SSZ = ZERO
176     !?              ENDIF
177     ! ----------------------------------------------------------------<<<
178     
179                 ELSE
180     
181     ! CARTESIAN GRID CASE
182     ! ---------------------------------------------------------------->>>
183     ! Surface forces
184     
185     ! bulk viscosity term
186                   SBV = (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     ! shear stress terms
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   ! CUT_CELL
214                     BCV = 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   ! same as NSW
236                             CUT_TAU_VS = .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   ! same as FSW
242                             CUT_TAU_VS = .FALSE.
243                             NOC_VS     = .FALSE.
244                             UW_s = ZERO
245                             VW_s = ZERO
246                             WW_s = ZERO
247                          ELSE                              ! partial slip
248                             CUT_TAU_VS = .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     ! SSX:
264                     U_NODE_AT_NE = ((.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     ! SSY:
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     ! SSZ:
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  ! end if do_k
383     
384                   ENDIF  ! end if/else cut_cell
385     ! ----------------------------------------------------------------<<<
386                 ENDIF   ! end if/else cartesian grid
387     
388     ! Source terms from shear stress
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     ! Add the terms
400                 TAU_V_S(IJK,M) = SBV + SSX + SSY + SSZ + Source_diff
401     
402               ELSE
403                 TAU_V_S(IJK,M) = ZERO
404               ENDIF   ! end if/else .not.sip_at_n .and. epsa>dil_ep_s
405             ENDDO   ! end do ijk
406     
407           ENDDO
408     
409           call send_recv(tau_v_s,2)
410           RETURN
411           END SUBROUTINE CALC_TAU_V_S
412