MFIX  2016-1
tau_w_g.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: CALC_Tau_W_g C
4 ! Purpose: Cross terms in the gradient of stress in W_g momentum C
5 ! C
6 ! Author: M. Syamlal Date: 18-DEC-96 C
7 ! C
8 ! C
9 ! Comments: This routine calculates the components of the gas phase C
10 ! viscous stress tensor of the w-momentum equation that cannot be C
11 ! cast in the form: mu.grad(w). These components are stored in the C
12 ! passed variable, which is then used as a source of the w-momentum C
13 ! equation. C
14 ! C
15 ! The following w component viscous stress tensor terms are C
16 ! calculated here: C
17 ! > part of 1/x d/dz (tau_zz) xdxdydz => C
18 ! 1/x d/dz (lambda.trcD) xdxdydz=> C
19 ! delta (lambda.trcD)Ap |T-B C
20 ! > part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently C
21 ! part of (tau_xz/x + 1/x d/dx (x tau_xz) ) xdxdydz => C
22 ! 1/x d/dx(mu.du/dz) xdxdydz => C
23 ! delta (mu/x du/dz)Ayz |E-W C
24 ! > part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently C
25 ! part of (tau_xz/x + 1/x d/dx (x tau_xz) ) xdxdydz => C
26 ! 1/x d/dx(mu.du/dz) xdxdydz => C
27 ! delta (mu/x du/dz)Ayz |E-W C
28 ! > part of 1/x d/dz (tau_zz) xdxdydz => C
29 ! 1/x d/dz (mu/x dw/dz) xdxdydz => C
30 ! delta (mu/x dw/dz)Axy |T-B C
31 ! CYLINDRICAL TERMS: C
32 ! > part of 1/x d/dz (tau_zz) xdxdydz => C
33 ! 1/x d/dz (2.mu/x u) xdxdydz => C
34 ! delta (2.mu/x u)Axy |T-B C
35 ! > part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently C
36 ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz => C
37 ! 1/x (mu/x du/dz) xdxdydz => C
38 ! delta (1/x mu/x du/dz)Vp C
39 ! C
40 ! The following w-component CYLINDRICAL viscous stress tensor terms C
41 ! are not calculated here. They are handled in source_w_g: C
42 ! > part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently C
43 ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz => C
44 ! 1/x d/dx (x.mu.(-w/x)) xdxdydz => C
45 ! delta (mu/x.(-w))Ayz |E-W C
46 ! > part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently C
47 ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz => C
48 ! mu/x dw/dx xdxdydz => C
49 ! delta (mu/x.(dw/dx))Vp |p C
50 ! > part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently C
51 ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz => C
52 ! 1/x d/dx (x.mu.(-w/x)) xdxdydz => C
53 ! delta (mu/x.(-w/x))Vp |p C
54 ! C
55 ! To reconstitute the complete w-momentum gas phase viscous stress C
56 ! tensor requires including the terms calculated in source_w_g C
57 ! and the 'diffusional' components (i.e., those of the form C
58 ! mu.grad(w) C
59 ! C
60 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
61  SUBROUTINE calc_tau_w_g(lTAU_W_G, lctau_w_g)
62 
63 ! Modules
64 !---------------------------------------------------------------------//
65  USE param
66  USE param1
67  USE constant
68  USE physprop
69  USE fldvar
70  USE visc_g
71  USE run
72  USE toleranc
73  USE geometry
74  USE indices
75  USE is
76  USE sendrecv
77  USE compar
78  USE functions
79  USE cutcell
80  IMPLICIT NONE
81 
82 ! Dummy arguments
83 !---------------------------------------------------------------------//
84 ! TAU_W_g
85  DOUBLE PRECISION, INTENT(OUT) :: lTAU_w_g(dimension_3)
86 ! cTAU_W_g
87  DOUBLE PRECISION, INTENT(OUT) :: lcTAU_w_g(dimension_3)
88 
89 ! Local variables
90 !---------------------------------------------------------------------//
91 ! Indices
92  INTEGER :: I, J, K, IM, JM, KP
93  INTEGER :: IJK, IJKE, IJKW, IJKN, IJKS, IJKT
94  INTEGER :: IJKNT, IJKST, IJKTE, IJKTW
95  INTEGER :: IJKP, IMJK, IJMK, IJKM
96  INTEGER :: IMJKP, IJMKP
97 ! Average volume fraction
98  DOUBLE PRECISION :: EPGA
99 ! Average gradients
100  DOUBLE PRECISION :: duodz
101 ! Source terms (Surface)
102  DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
103 ! Source terms (Volumetric)
104  DOUBLE PRECISION :: Vxz
105 !---------------------------------------------------------------------//
106 
107  IF((.NOT.cartesian_grid).OR.(cg_safe_mode(5)==1)) THEN
108 
109 !$omp parallel do default(none) &
110 !$omp private(I, J, K, IM, JM, KP, IJK, &
111 !$omp IJKE, IJKW, IJKN, IJKS, IJKT, &
112 !$omp IJKTE, IJKNT, IJKTW, IJKST, &
113 !$omp IMJK, IJMK, IJKM, IJKP, IJMKP, IMJKP, &
114 !$omp EPGA, SBV, SSX, SSY, SSZ, vxz, duodz) &
115 !$omp shared(ijkstart3, ijkend3, i_of, j_of, k_of, im1, jm1, kp1, &
116 !$omp do_k, cylindrical, ltau_w_g, lctau_w_g, &
117 !$omp ep_g, epmu_gt, lambda_gt, trd_g, v_g, w_g, u_g, &
118 !$omp axy, axy_w, ayz_w, axz_w, vol_w, &
119 !$omp dy, dz, ox, ox_e, odz, odz_t, eplambda_gt)
120  DO ijk = ijkstart3, ijkend3
121  k = k_of(ijk)
122  ijkt = top_of(ijk)
123  epga = avg_z(ep_g(ijk),ep_g(ijkt),k)
124  IF ( .NOT.ip_at_t(ijk) .AND. epga>dil_ep_s) THEN
125  j = j_of(ijk)
126  i = i_of(ijk)
127  im = im1(i)
128  jm = jm1(j)
129  kp = kp1(k)
130 
131  ijkp = kp_of(ijk)
132  imjk = im_of(ijk)
133  ijmk = jm_of(ijk)
134  ijkm = km_of(ijk)
135  imjkp = kp_of(imjk)
136  ijmkp = jm_of(ijkp)
137 
138  ijkn = north_of(ijk)
139  ijks = south_of(ijk)
140  ijke = east_of(ijk)
141  ijkw = west_of(ijk)
142  ijknt = top_of(ijkn)
143  ijkst = top_of(ijks)
144  ijkte = east_of(ijkt)
145  ijktw = west_of(ijkt)
146 
147 
148 ! Surface forces
149 ! Bulk viscosity term
150 ! part of 1/x d/dz (tau_zz) xdxdydz =>
151 ! 1/x d/dz (lambda.trcD) xdxdydz=>
152 ! delta (lambda.trcD)Ap |T-B : at (i, j, k+1 - k-1)
153  sbv = (eplambda_gt(ijkt)*trd_g(ijkt)-&
154  eplambda_gt(ijk)*trd_g(ijk))*axy(ijk)
155 
156 ! shear stress terms
157 ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
158 ! part of (tau_xz/x + 1/x d/dx (x tau_xz) ) xdxdydz =>
159 ! 1/x d/dx(mu.du/dz) xdxdydz =>
160 ! delta (mu/x du/dz)Ayz |E-W : at (i+1/2-i-1/2, j, k+1/2)
161  ssx = avg_z_h(avg_x_h(epmu_gt(ijk),epmu_gt(ijke),i),&
162  avg_x_h(epmu_gt(ijkt),epmu_gt(ijkte),i),k)*&
163  (u_g(ijkp)-u_g(ijk))*ox_e(i)*odz_t(k)*&
164  ayz_w(ijk) - &
165  avg_z_h(avg_x_h(epmu_gt(ijkw),epmu_gt(ijk),im),&
166  avg_x_h(epmu_gt(ijktw),epmu_gt(ijkt),im),k)*&
167  (u_g(imjkp)-u_g(imjk))*odz_t(k)*dy(j)*&
168  (half*(dz(k)+dz(kp)))
169 ! DY(J)*HALF(DZ(k)+DZ(kp)) = oX_E(IM)*AYZ_W(IMJK), but avoids singularity
170 
171 ! part of d/dy (tau_zy) xdxdydz =>
172 ! d/dy (mu/x dv/dz) xdxdydz =>
173 ! delta (mu/x dv/dz)Axz |N-S : at (i, j+1/2 - j-1/2, k+1/2)
174  ssy = avg_z_h(avg_y_h(epmu_gt(ijk),epmu_gt(ijkn),j),&
175  avg_y_h(epmu_gt(ijkt),epmu_gt(ijknt),j),k)*&
176  (v_g(ijkp)-v_g(ijk))*ox(i)*odz_t(k)*axz_w(ijk) -&
177  avg_z_h(avg_y_h(epmu_gt(ijks),epmu_gt(ijk),jm),&
178  avg_y_h(epmu_gt(ijkst),epmu_gt(ijkt),jm),k)*&
179  (v_g(ijmkp)-v_g(ijmk))*ox(i)*odz_t(k)*axz_w(ijmk)
180 
181 ! part of 1/x d/dz (tau_zz) xdxdydz =>
182 ! 1/x d/dz (mu/x dw/dz) xdxdydz =>
183 ! delta (mu/x dw/dz)Axy |T-B : at (i, j, k+1 - k-1)
184  ssz = epmu_gt(ijkt)*(w_g(ijkp)-w_g(ijk))*ox(i)*odz(kp)*&
185  axy_w(ijk) - &
186  epmu_gt(ijk)*(w_g(ijk)-w_g(ijkm))*ox(i)*odz(k)*&
187  axy_w(ijkm)
188 
189 
190 ! Special terms for cylindrical coordinates
191  IF (cylindrical) THEN
192 
193 ! part of 1/x d/dz (tau_zz) xdxdydz =>
194 ! 1/x d/dz (2.mu/x u) xdxdydz =>
195 ! delta (2.mu/x u)Axy |T-B : at (i, j, k+1 - k-1)
196  ssz = ssz + epmu_gt(ijkt)*(u_g(ijkp)+u_g(imjkp))*&
197  ox(i)*axy_w(ijk) - &
198  epmu_gt(ijk)*(u_g(ijk)+u_g(imjk))*&
199  ox(i)*axy_w(ijkm)
200 
201 ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
202 ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
203 ! 1/x (mu/x du/dz) xdxdydz =>
204 ! delta (1/x mu/x du/dz)Vp : at (i, j, k+1/2)
205  IF (ox_e(im) /= undefined) THEN
206  duodz = (u_g(imjkp)-u_g(imjk))*ox_e(im)*odz_t(k)
207  ELSE
208  duodz = zero
209  ENDIF
210  vxz = avg_z(epmu_gt(ijk),epmu_gt(ijkt),k)*ox(i)*half*&
211  ( (u_g(ijkp)-u_g(ijk))*ox_e(i)*odz_t(k) + &
212  duodz )
213  ELSE
214  vxz = zero
215  ENDIF
216 
217 ! Add the terms
218  ltau_w_g(ijk) = sbv + ssx + ssy + ssz + vxz*vol_w(ijk)
219 
220 ! Also calculate and store the full gas phase viscous stress tensor
221  CALL get_full_tau_w_g(ijk, ltau_w_g, lctau_w_g)
222 
223  ELSE
224  ltau_w_g(ijk) = zero
225  lctau_w_g(ijk) = zero
226  ENDIF ! end if (.NOT. IP_AT_T(IJK) .AND. EPGA>DIL_EP_S)
227  ENDDO ! end do ijk
228 !$omp end parallel do
229 
230  ELSE
231 ! if cartesian grid
232  CALL calc_cg_tau_w_g(ltau_w_g, lctau_w_g)
233  ENDIF
234 
235  call send_recv(ltau_w_g,2)
236  call send_recv(lctau_w_g,2)
237 
238  RETURN
239 
240  CONTAINS
241 
242  include 'fun_avg.inc'
243 
244  END SUBROUTINE calc_tau_w_g
245 
246 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
247 ! C
248 ! Subroutine: CALC_CG_Tau_W_g C
249 ! Purpose: Cross terms in the gradient of stress in w_g momentum C
250 ! based on cartesian grid cut cell. C
251 ! C
252 ! Author: Jeff Dietiker Date: 01-Jul-09 C
253 ! C
254 ! C
255 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
256  SUBROUTINE calc_cg_tau_w_g(lTAu_w_G, lctau_w_g)
257 
258 ! Modules
259 !---------------------------------------------------------------------//
260  USE param
261  USE param1
262  USE constant
263  USE physprop
264  USE fldvar
265  USE visc_g
266  USE run
267  USE toleranc
268  USE geometry
269  USE indices
270  USE is
271  USE sendrecv
272  USE compar
273  USE fun_avg
274  USE functions
275 
276  USE bc
277  USE quadric
278  USE cutcell
279  IMPLICIT NONE
280 
281 ! Dummy arguments
282 !---------------------------------------------------------------------//
283 ! TAU_W_g
284  DOUBLE PRECISION, INTENT(OUT) :: lTAU_w_g(dimension_3)
285 ! cTAU_W_g
286  DOUBLE PRECISION, INTENT(OUT) :: lcTAU_w_g(dimension_3)
287 
288 ! Local variables
289 !---------------------------------------------------------------------//
290 ! Indices
291  INTEGER :: I, J, K, IM, JM, KP
292  INTEGER :: IJK, IJKE, IJKW, IJKN, IJKS, IJKT
293  INTEGER :: IJKNT, IJKST, IJKTE, IJKTW
294  INTEGER :: IJKP, IMJK, IJMK, IJKM
295  INTEGER :: IMJKP, IJMKP
296 ! Average volume fraction
297  DOUBLE PRECISION :: EPGA
298 ! Source terms (Surface)
299  DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
300 ! cartesian grid
301  DOUBLE PRECISION :: DEL_H, Nx, Ny, Nz
302  LOGICAL :: U_NODE_AT_ET, U_NODE_AT_EB, U_NODE_AT_WT, U_NODE_AT_WB
303  LOGICAL :: V_NODE_AT_NT, V_NODE_AT_NB, V_NODE_AT_ST, V_NODE_AT_SB
304  DOUBLE PRECISION :: dudz_at_E, dudz_at_W
305  DOUBLE PRECISION :: dvdz_at_N, dvdz_at_S
306  DOUBLE PRECISION :: Xi, Yi, Zi, Ui, Vi, Sx, Sy, Sz
307  DOUBLE PRECISION :: MU_GT_CUT, SSX_CUT, SSY_CUT
308  DOUBLE PRECISION :: UW_g, VW_g, WW_g
309  INTEGER :: BCV
310  INTEGER :: BCT
311 
312 !---------------------------------------------------------------------//
313 
314 !$omp parallel do default(none) &
315 !$omp private(I, J, K, IM, JM, KP, IJK, &
316 !$omp IJKE, IJKW, IJKN, IJKS, IJKT, &
317 !$omp IJKTE, IJKNT, IJKTW, IJKST, &
318 !$omp IMJK, IJMK, IJKM, IJKP, IJMKP, IMJKP, &
319 !$omp EPGA, SBV, SSX, SSY, SSZ, &
320 !$omp BCV, BCT, BC_TYPE_ENUM, noc_wg, uw_g, vw_g, ww_g, &
321 !$omp del_h, nx, ny, nz, xi, yi, zi, sx, sy, sz, vi, ui, &
322 !$omp cut_tau_wg, mu_gt_cut, ssy_cut, ssx_cut, &
323 !$omp dudz_at_e, dudz_at_w, dvdz_at_S, dvdz_at_N, &
324 !$omp u_node_at_wb, u_node_at_wt, u_node_at_eb, u_node_at_et, &
325 !$omp v_node_at_sb, v_node_at_st, v_node_at_nb, v_node_at_nt) &
326 !$omp shared(ijkstart3, ijkend3, i_of, j_of, k_of, im1, jm1, kp1, &
327 !$omp do_k, cylindrical, ltau_w_g, lctau_w_g, &
328 !$omp ep_g, mu_g, epmu_gt, lambda_gt, trd_g, v_g, w_g, u_g, &
329 !$omp axy, axy_w, ayz_w, axz_w, vol, ox, &
330 !$omp bc_type, bc_w_id, bc_hw_g, bc_uw_g, bc_vw_g, bc_ww_g, &
331 !$omp oneodz_t_u, oneodz_t_v, oneodz_t_w, &
332 !$omp x_u, y_u, z_u, x_v, y_v, z_v, x_w, y_w, z_w, &
333 !$omp wall_u_at, wall_v_at, area_w_cut, cut_w_cell_at, &
334 !$omp blocked_u_cell_at, blocked_v_cell_at, eplambda_gt)
335  DO ijk = ijkstart3, ijkend3
336  k = k_of(ijk)
337  ijkt = top_of(ijk)
338  epga = avg_z(ep_g(ijk),ep_g(ijkt),k)
339  IF ( .NOT.ip_at_t(ijk) .AND. epga>dil_ep_s) THEN
340  j = j_of(ijk)
341  i = i_of(ijk)
342  im = im1(i)
343  jm = jm1(j)
344  kp = kp1(k)
345 
346  ijkp = kp_of(ijk)
347  imjk = im_of(ijk)
348  ijmk = jm_of(ijk)
349  ijkm = km_of(ijk)
350  imjkp = kp_of(imjk)
351  ijmkp = jm_of(ijkp)
352 
353  ijkn = north_of(ijk)
354  ijks = south_of(ijk)
355  ijke = east_of(ijk)
356  ijkw = west_of(ijk)
357  ijknt = top_of(ijkn)
358  ijkst = top_of(ijks)
359  ijkte = east_of(ijkt)
360  ijktw = west_of(ijkt)
361 
362 ! bulk viscosity term
363  sbv = (eplambda_gt(ijkt)*trd_g(ijkt)) * axy_w(ijk) &
364  -(eplambda_gt(ijk) *trd_g(ijk) ) * axy_w(ijkm)
365 
366 ! shear stress terms
367  IF(.NOT.cut_w_cell_at(ijk)) THEN
368  ssx = avg_z_h(avg_x_h(epmu_gt(ijk),epmu_gt(ijke),i),&
369  avg_x_h(epmu_gt(ijkt),epmu_gt(ijkte),i),k)*&
370  (u_g(ijkp)-u_g(ijk))*oneodz_t_u(ijk)*ayz_w(ijk) -&
371  avg_z_h(avg_x_h(epmu_gt(ijkw),epmu_gt(ijk),im),&
372  avg_x_h(epmu_gt(ijktw),epmu_gt(ijkt),im),k)*&
373  (u_g(imjkp)-u_g(imjk))*oneodz_t_u(imjk)*ayz_w(imjk)
374 
375  ssy = avg_z_h(avg_y_h(epmu_gt(ijk),epmu_gt(ijkn),j),&
376  avg_y_h(epmu_gt(ijkt),epmu_gt(ijknt),j),k)*&
377  (v_g(ijkp)-v_g(ijk))*oneodz_t_v(ijk)*axz_w(ijk) - &
378  avg_z_h(avg_y_h(epmu_gt(ijks),epmu_gt(ijk),jm),&
379  avg_y_h(epmu_gt(ijkst),epmu_gt(ijkt),jm),k)*&
380  (v_g(ijmkp)-v_g(ijmk))*oneodz_t_v(ijmk)*axz_w(ijmk)
381 
382  ssz = epmu_gt(ijkt)*(w_g(ijkp)-w_g(ijk))*&
383  oneodz_t_w(ijk)*axy_w(ijk) - &
384  epmu_gt(ijk)*(w_g(ijk)-w_g(ijkm))*&
385  oneodz_t_w(ijkm)*axy_w(ijkm)
386 
387 ! cut cell modifications
388 !---------------------------------------------------------------------//
389  ELSE
390  bcv = bc_w_id(ijk)
391 
392  IF(bcv > 0 ) THEN
393  bct = bc_type_enum(bcv)
394  ELSE
395  bct = none
396  ENDIF
397 
398  SELECT CASE (bct)
399  CASE (cg_nsw,cg_mi)
400  cut_tau_wg = .true.
401  noc_wg = .true.
402  uw_g = zero
403  vw_g = zero
404  ww_g = zero
405  CASE (cg_fsw)
406  cut_tau_wg = .false.
407  noc_wg = .false.
408  uw_g = zero
409  vw_g = zero
410  ww_g = zero
411  CASE(cg_psw)
412  IF(bc_hw_g(bc_w_id(ijk))==undefined) THEN ! same as NSW
413  cut_tau_wg = .true.
414  noc_wg = .true.
415  uw_g = bc_uw_g(bcv)
416  vw_g = bc_vw_g(bcv)
417  ww_g = bc_ww_g(bcv)
418  ELSEIF(bc_hw_g(bc_w_id(ijk))==zero) THEN ! same as FSW
419  cut_tau_wg = .false.
420  noc_wg = .false.
421  uw_g = zero
422  vw_g = zero
423  ww_g = zero
424  ELSE ! partial slip
425  cut_tau_wg = .false.
426  noc_wg = .false.
427  ENDIF
428  CASE (none)
429  ltau_w_g(ijk) = zero
430  cycle
431  END SELECT
432 
433  IF(cut_tau_wg) THEN
434  mu_gt_cut = (vol(ijk)*epmu_gt(ijk) + &
435  vol(ijkp)*epmu_gt(ijkt))/(vol(ijk) + vol(ijkp))
436  ELSE
437  mu_gt_cut = zero
438  ENDIF
439 
440 ! SSX:
441  u_node_at_et = ((.NOT.blocked_u_cell_at(ijkp)).AND.&
442  (.NOT.wall_u_at(ijkp)))
443  u_node_at_eb = ((.NOT.blocked_u_cell_at(ijk)).AND.&
444  (.NOT.wall_u_at(ijk)))
445  u_node_at_wt = ((.NOT.blocked_u_cell_at(imjkp)).AND.&
446  (.NOT.wall_u_at(imjkp)))
447  u_node_at_wb = ((.NOT.blocked_u_cell_at(imjk)).AND.&
448  (.NOT.wall_u_at(imjk)))
449 
450  IF(u_node_at_et.AND.u_node_at_eb) THEN
451  ui = half * (u_g(ijkp) + u_g(ijk))
452  xi = half * (x_u(ijkp) + x_u(ijk))
453  yi = half * (y_u(ijkp) + y_u(ijk))
454  zi = half * (z_u(ijkp) + z_u(ijk))
455  sx = x_u(ijkp) - x_u(ijk)
456  sy = y_u(ijkp) - y_u(ijk)
457  sz = z_u(ijkp) - z_u(ijk)
458  CALL get_del_h(ijk, 'W_MOMENTUM', xi, yi, zi, &
459  del_h, nx, ny, nz)
460 
461  dudz_at_e = (u_g(ijkp) - u_g(ijk)) * &
462  oneodz_t_u(ijk)
463  IF(noc_wg) dudz_at_e = dudz_at_e - ((ui-uw_g) * &
464  oneodz_t_u(ijk)/del_h*(sx*nx+sy*ny))
465  ELSE
466  dudz_at_e = zero
467  ENDIF
468 
469  IF(u_node_at_wt.AND.u_node_at_wb) THEN
470  ui = half * (u_g(imjkp) + u_g(imjk))
471  xi = half * (x_u(imjkp) + x_u(imjk))
472  yi = half * (y_u(imjkp) + y_u(imjk))
473  zi = half * (z_u(imjkp) + z_u(imjk))
474  sx = x_u(imjkp) - x_u(imjk)
475  sy = y_u(imjkp) - y_u(imjk)
476  sz = z_u(imjkp) - z_u(imjk)
477  CALL get_del_h(ijk, 'W_MOMENTUM', xi, yi, zi, &
478  del_h, nx, ny, nz)
479 
480  dudz_at_w = (u_g(imjkp) - u_g(imjk)) * &
481  oneodz_t_u(imjk)
482  IF(noc_wg) dudz_at_w = dudz_at_w - ((ui-uw_g) * &
483  oneodz_t_u(imjk)/del_h*(sx*nx+sy*ny))
484  ELSE
485  dudz_at_w = zero
486  ENDIF
487 
488  IF(u_node_at_eb) THEN
489  CALL get_del_h(ijk, 'W_MOMENTUM', x_u(ijk), &
490  y_u(ijk), z_u(ijk), del_h, nx, ny, nz)
491  ssx_cut = - mu_gt_cut * (u_g(ijk) - uw_g) / &
492  del_h * (nz*nx) * area_w_cut(ijk)
493  ELSE
494  ssx_cut = zero
495  ENDIF
496  ssx = avg_z_h(avg_x_h(epmu_gt(ijk),epmu_gt(ijke),i),&
497  avg_x_h(epmu_gt(ijkt),epmu_gt(ijkte),i),k)*&
498  dudz_at_e*ayz_w(ijk) - &
499  avg_z_h(avg_x_h(epmu_gt(ijkw),epmu_gt(ijk),im),&
500  avg_x_h(epmu_gt(ijktw),epmu_gt(ijkt),im),k)*&
501  dudz_at_w*ayz_w(imjk) + ssx_cut
502 
503 ! SSY:
504  v_node_at_nt = ((.NOT.blocked_v_cell_at(ijkp)).AND.&
505  (.NOT.wall_v_at(ijkp)))
506  v_node_at_nb = ((.NOT.blocked_v_cell_at(ijk)).AND.&
507  (.NOT.wall_v_at(ijk)))
508  v_node_at_st = ((.NOT.blocked_v_cell_at(ijmkp)).AND.&
509  (.NOT.wall_v_at(ijmkp)))
510  v_node_at_sb = ((.NOT.blocked_v_cell_at(ijmk)).AND.&
511  (.NOT.wall_v_at(ijmk)))
512 
513  IF(v_node_at_nt.AND.v_node_at_nb) THEN
514  vi = half * (v_g(ijkp) + v_g(ijk))
515  xi = half * (x_v(ijkp) + x_v(ijk))
516  yi = half * (y_v(ijkp) + y_v(ijk))
517  zi = half * (z_v(ijkp) + z_v(ijk))
518  sx = x_v(ijkp) - x_v(ijk)
519  sy = y_v(ijkp) - y_v(ijk)
520  sz = z_v(ijkp) - z_v(ijk)
521  CALL get_del_h(ijk, 'W_MOMENTUM', xi, yi, zi, &
522  del_h, nx, ny, nz)
523 
524  dvdz_at_n = (v_g(ijkp) - v_g(ijk)) * &
525  oneodz_t_v(ijk)
526  IF(noc_wg) dvdz_at_n = dvdz_at_n - ((vi-vw_g) * &
527  oneodz_t_v(ijk)/del_h*(sx*nx+sy*ny))
528  ELSE
529  dvdz_at_n = zero
530  ENDIF
531 
532  IF(v_node_at_st.AND.v_node_at_sb) THEN
533  vi = half * (v_g(ijmkp) + v_g(ijmk))
534  xi = half * (x_v(ijmkp) + x_v(ijmk))
535  yi = half * (y_v(ijmkp) + y_v(ijmk))
536  zi = half * (z_v(ijmkp) + z_v(ijmk))
537  sx = x_v(ijmkp) - x_v(ijmk)
538  sy = y_v(ijmkp) - y_v(ijmk)
539  sz = z_v(ijmkp) - z_v(ijmk)
540  CALL get_del_h(ijk, 'W_MOMENTUM', xi, yi, zi, &
541  del_h, nx, ny, nz)
542 
543  dvdz_at_s = (v_g(ijmkp) - v_g(ijmk)) * &
544  oneodz_t_v(ijmk)
545  IF(noc_wg) dvdz_at_s = dvdz_at_s - ((vi-vw_g) * &
546  oneodz_t_v(ijmk)/del_h*(sx*nx+sy*ny))
547  ELSE
548  dvdz_at_s = zero
549  ENDIF
550 
551  IF(v_node_at_nb) THEN
552  CALL get_del_h(ijk, 'W_MOMENTUM', x_v(ijk), &
553  y_v(ijk), z_v(ijk), del_h, nx, ny, nz)
554  ssy_cut = - mu_gt_cut * (v_g(ijk) - vw_g) / &
555  del_h * (nz*ny) * area_w_cut(ijk)
556  ELSE
557  ssy_cut = zero
558  ENDIF
559 
560  ssy = avg_z_h(avg_y_h(epmu_gt(ijk),epmu_gt(ijkn),j),&
561  avg_y_h(epmu_gt(ijkt),epmu_gt(ijknt),j),k)*&
562  dvdz_at_n*axz_w(ijk) - &
563  avg_z_h(avg_y_h(epmu_gt(ijks),epmu_gt(ijk),jm),&
564  avg_y_h(epmu_gt(ijkst),epmu_gt(ijkt),jm),k)*&
565  dvdz_at_s*axz_w(ijmk) + ssy_cut
566 
567 ! SSZ:
568  CALL get_del_h(ijk, 'W_MOMENTUM', x_w(ijk), &
569  y_w(ijk), z_w(ijk), del_h, nx, ny, nz)
570 
571  ssz = epmu_gt(ijkt)*(w_g(ijkp)-w_g(ijk))*&
572  oneodz_t_w(ijk)*axy_w(ijk) - &
573  epmu_gt(ijk)*(w_g(ijk)-w_g(ijkm))*&
574  ox(i)*oneodz_t_w(ijkm)*axy_w(ijkm) - &
575  mu_gt_cut * (w_g(ijk) - ww_g) / del_h * &
576  (nz**2) * area_w_cut(ijk)
577 
578  ENDIF ! end if/else cut_w_cell_at
579 
580 ! Add the terms
581  ltau_w_g(ijk) = sbv + ssx + ssy + ssz
582 
583 ! Also calculate and store the full gas phase viscous stress tensor
584  CALL get_full_tau_w_g(ijk, ltau_w_g, lctau_w_g)
585 
586  ELSE
587  ltau_w_g(ijk) = zero
588  lctau_w_g(ijk) = zero
589  ENDIF ! end if (.NOT. IP_AT_T(IJK) .AND. EPGA>DIL_EP_S)
590  ENDDO ! end do ijk
591 !$omp end parallel do
592 
593  RETURN
594  END SUBROUTINE calc_cg_tau_w_g
595 
596 
597 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
598 ! C
599 ! Subroutine: GET_FULL_Tau_W_g C
600 ! Purpose: Calculate the divergence of the complete gas phase stress C
601 ! tensor including all terms in the w-direction. C
602 ! C
603 ! C
604 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
605  SUBROUTINE get_full_tau_w_g(IJK, ltau_w_g, lctau_w_g)
606 
607 ! Modules
608 !---------------------------------------------------------------------//
609  USE fldvar, only: w_g
610 
611  USE functions, only: east_of, west_of, top_of, flow_at_t
612  USE functions, only: im_of, jm_of, km_of
613  USE functions, only: ip_of, jp_of, kp_of
614 
615  USE fun_avg, only: avg_x_h, avg_z_h, avg_z
616 
617  USE geometry, only: cylindrical, ox_e, ox, odx_e
618  USE geometry, only: dy, dz, vol_w, ayz_w
619  USE indices, only: i_of, j_of, k_of, im1, kp1
620 
621  use param
622  USE param1, only: zero, half
623  USE visc_g, only: mu_gt, df_gw
624  IMPLICIT NONE
625 
626 ! Dummy arguments
627 !---------------------------------------------------------------------//
628 ! ijk index
629  INTEGER, INTENT(IN) :: IJK
630 ! TAU_W_g
631  DOUBLE PRECISION, INTENT(IN) :: lTAU_W_g(dimension_3)
632 ! cTAU_W_g
633  DOUBLE PRECISION, INTENT(INOUT) :: lcTAU_W_g(dimension_3)
634 
635 ! Local Variables
636 !---------------------------------------------------------------------//
637 ! indices
638  INTEGER :: I, J, K, IM
639  INTEGER :: IJKT, IJKE, IJKW, IJKTE, IJKTW
640  INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
641 ! average viscosity
642  DOUBLE PRECISION :: MUOX
643 ! source terms
644  DOUBLE PRECISION :: SSX, SSY, SSZ
645 ! cylindrical terms
646  DOUBLE PRECISION :: vxza
647  DOUBLE PRECISION :: cte, ctw, vx14
648  DOUBLE PRECISION :: cpe, cpw, vx15
649 !---------------------------------------------------------------------//
650 
651  ipjk = ip_of(ijk)
652  ijpk = jp_of(ijk)
653  ijkp = kp_of(ijk)
654  imjk = im_of(ijk)
655  ijmk = jm_of(ijk)
656  ijkm = km_of(ijk)
657 
658 ! additional terms (found in source_W_g)
659  vxza = zero
660  vx14 = zero
661  vx15 = zero
662  IF (cylindrical) THEN
663  i = i_of(ijk)
664  j = j_of(ijk)
665  k = k_of(ijk)
666  im = im1(i)
667  ijkt = top_of(ijk)
668  ijke = east_of(ijk)
669  ijkw = west_of(ijk)
670  ijkte = top_of(ijke)
671  ijktw = top_of(ijkw)
672 
673 ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
674 ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
675 ! 1/x d/dx (x.mu.(-w/x)) xdxdydz =>
676 ! delta (mu/x.(-w))Ayz |E-W : at (i+1/2 - i-1/2, j, k+1/2)
677  cte = half*avg_z_h(avg_x_h(mu_gt(ijk),mu_gt(ijke),i),&
678  avg_x_h(mu_gt(ijkt),mu_gt(ijkte),i),k)*&
679  ox_e(i)*ayz_w(ijk)
680  ctw = half*avg_z_h(avg_x_h(mu_gt(ijkw),mu_gt(ijk),im),&
681  avg_x_h(mu_gt(ijktw),mu_gt(ijkt),im),k)*&
682  dy(j)*(half*(dz(k)+dz(kp1(k))))
683 ! DY(J)*HALF(DZ(k)+DZ(kp)) = oX_E(IM)*AYZ_W(IMJK), but avoids singularity
684  vx14 = -cte*(w_g(ipjk)+w_g(ijk))+ctw*(w_g(ijk)+w_g(imjk))
685 
686 ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
687 ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
688 ! mu/x dw/dx xdxdydz =>
689 ! delta (mu/x.(dw/dx))Vp |p : at (i, j, k+1/2)
690  muox = avg_z(mu_gt(ijk),mu_gt(ijkt),k)*ox(i)
691  cpe = half*muox*odx_e(i)*vol_w(ijk)
692  cpw = half*muox*odx_e(im)*vol_w(ijk)
693  vx15 = cpe*(w_g(ipjk)-w_g(ijk))+cpw*(w_g(ijk)-w_g(imjk))
694 
695 ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
696 ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
697 ! 1/x d/dx (x.mu.(-w/x)) xdxdydz =>
698 ! delta (mu/x.(-w/x))Vp |p : at (i, j, k+1/2)
699  vxza = -muox*ox(i)*w_g(ijk)*vol_w(ijk)
700  ENDIF
701 
702 
703 ! convection terms
704  ssx = zero
705  ssy = zero
706  ssz = zero
707  IF (flow_at_t(ijk)) THEN
708 ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
709 ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
710 ! 1/x d/dx (x.mu.dw/dx) xdxdydz =>
711 ! delta (mu.dw/dx.Ayz) |E-W : at (i+1/2 - i-1/2), j, k+1/2
712  ssx = df_gw(ijk,east)*(w_g(ipjk) - w_g(ijk)) - &
713  df_gw(ijk,west)*(w_g(ijk) - w_g(ijkm))
714 
715 ! part of d/dy (tau_zy) xdxdydz =>
716 ! d/dy (mu.dw/dy) xdxdydz =>
717 ! delta (mu.dw/dy.Axz) |N-S : at (i, j+1/2 - j-1/2, k+1/2)
718  ssy = df_gw(ijk,north)*(w_g(ijpk)-w_g(ijk)) - &
719  df_gw(ijk,south)*(w_g(ijk)-w_g(ijmk))
720 
721 ! part of 1/x d/dz (tau_zz) xdxdydz =>
722 ! 1/x d/dz (mu/x.dw/dz) xdxdydz =>
723 ! delta (mu/x.dw/dz.Axy) |T-B : at (i, j, k+1 - k-1)
724  ssz = df_gw(ijk,top)*(w_g(ijkp)-w_g(ijk)) - &
725  df_gw(ijk,bottom)*(w_g(ijk)-w_g(ijkm))
726 
727  ENDIF
728 
729 ! Add the terms
730  lctau_w_g(ijk) = (ltau_w_g(ijk) + ssx + ssy + ssz + vxza + &
731  vx14 + vx15)
732 
733  RETURN
734  END SUBROUTINE get_full_tau_w_g
735 
cutcell::cut_tau_wg
logical cut_tau_wg
Definition: cutcell_mod.f:390
cutcell::y_v
double precision, dimension(:), allocatable y_v
Definition: cutcell_mod.f:55
geometry::vol_w
double precision, dimension(:), allocatable vol_w
Definition: geometry_mod.f:242
cutcell::z_u
double precision, dimension(:), allocatable z_u
Definition: cutcell_mod.f:51
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
cutcell::wall_u_at
logical, dimension(:), allocatable wall_u_at
Definition: cutcell_mod.f:126
bc::bc_uw_g
double precision, dimension(dimension_bc) bc_uw_g
Definition: bc_mod.f:313
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
fldvar::ep_g
double precision, dimension(:), allocatable ep_g
Definition: fldvar_mod.f:15
visc_g::df_gw
double precision, dimension(:,:), allocatable df_gw
Definition: visc_g_mod.f:33
geometry::ox_e
double precision, dimension(:), allocatable ox_e
Definition: geometry_mod.f:143
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
visc_g::mu_gt
double precision, dimension(:), allocatable mu_gt
Definition: visc_g_mod.f:8
geometry::axy
double precision, dimension(:), allocatable axy
Definition: geometry_mod.f:210
cutcell::oneodz_t_w
double precision, dimension(:), allocatable oneodz_t_w
Definition: cutcell_mod.f:326
cutcell::x_u
double precision, dimension(:), allocatable x_u
Definition: cutcell_mod.f:49
indices::im1
integer, dimension(:), allocatable im1
Definition: indices_mod.f:50
cutcell::cg_safe_mode
integer, dimension(10) cg_safe_mode
Definition: cutcell_mod.f:415
cutcell::wall_v_at
logical, dimension(:), allocatable wall_v_at
Definition: cutcell_mod.f:127
geometry::dy
double precision, dimension(0:dim_j) dy
Definition: geometry_mod.f:70
bc::bc_type_enum
integer, dimension(dimension_bc) bc_type_enum
Definition: bc_mod.f:146
param1::undefined
double precision, parameter undefined
Definition: param1_mod.f:18
param::east
integer east
Definition: param_mod.f:29
cutcell::y_w
double precision, dimension(:), allocatable y_w
Definition: cutcell_mod.f:60
geometry::dz
double precision, dimension(0:dim_k) dz
Definition: geometry_mod.f:72
cutcell::z_v
double precision, dimension(:), allocatable z_v
Definition: cutcell_mod.f:56
calc_cg_tau_w_g
subroutine calc_cg_tau_w_g(lTAu_w_G, lctau_w_g)
Definition: tau_w_g.f:257
is
Definition: is_mod.f:11
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
fun_avg
Definition: fun_avg_mod.f:1
sendrecv::send_recv
Definition: sendrecv_mod.f:75
cutcell::noc_wg
logical noc_wg
Definition: cutcell_mod.f:389
cutcell::x_w
double precision, dimension(:), allocatable x_w
Definition: cutcell_mod.f:59
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
geometry::odx_e
double precision, dimension(:), allocatable odx_e
Definition: geometry_mod.f:121
indices::jm1
integer, dimension(:), allocatable jm1
Definition: indices_mod.f:51
visc_g::trd_g
double precision, dimension(:), allocatable trd_g
Definition: visc_g_mod.f:4
cutcell::bc_w_id
integer, dimension(:), allocatable bc_w_id
Definition: cutcell_mod.f:436
geometry::ox
double precision, dimension(:), allocatable ox
Definition: geometry_mod.f:140
visc_g::epmu_gt
double precision, dimension(:), allocatable epmu_gt
Definition: visc_g_mod.f:13
cutcell::blocked_u_cell_at
logical, dimension(:), allocatable blocked_u_cell_at
Definition: cutcell_mod.f:364
bc::bc_hw_g
double precision, dimension(dimension_bc) bc_hw_g
Definition: bc_mod.f:307
param::north
integer north
Definition: param_mod.f:37
fldvar::v_g
double precision, dimension(:), allocatable v_g
Definition: fldvar_mod.f:99
cutcell::cut_w_cell_at
logical, dimension(:), allocatable cut_w_cell_at
Definition: cutcell_mod.f:358
param::south
integer south
Definition: param_mod.f:41
indices::kp1
integer, dimension(:), allocatable kp1
Definition: indices_mod.f:52
fldvar::w_g
double precision, dimension(:), allocatable w_g
Definition: fldvar_mod.f:111
cutcell::x_v
double precision, dimension(:), allocatable x_v
Definition: cutcell_mod.f:54
param1::half
double precision, parameter half
Definition: param1_mod.f:28
run
Definition: run_mod.f:13
geometry::ayz_w
double precision, dimension(:), allocatable ayz_w
Definition: geometry_mod.f:236
cutcell::oneodz_t_u
double precision, dimension(:), allocatable oneodz_t_u
Definition: cutcell_mod.f:317
param
Definition: param_mod.f:2
toleranc::dil_ep_s
double precision, parameter dil_ep_s
Definition: toleranc_mod.f:24
cutcell::cartesian_grid
logical cartesian_grid
Definition: cutcell_mod.f:13
visc_g::eplambda_gt
double precision, dimension(:), allocatable eplambda_gt
Definition: visc_g_mod.f:21
geometry::axz_w
double precision, dimension(:), allocatable axz_w
Definition: geometry_mod.f:238
geometry::odz
double precision, dimension(:), allocatable odz
Definition: geometry_mod.f:118
geometry::cylindrical
logical cylindrical
Definition: geometry_mod.f:168
bc::bc_vw_g
double precision, dimension(dimension_bc) bc_vw_g
Definition: bc_mod.f:316
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
visc_g
Definition: visc_g_mod.f:1
param::west
integer west
Definition: param_mod.f:33
geometry::axy_w
double precision, dimension(:), allocatable axy_w
Definition: geometry_mod.f:240
get_full_tau_w_g
subroutine get_full_tau_w_g(IJK, ltau_w_g, lctau_w_g)
Definition: tau_w_g.f:606
fldvar::u_g
double precision, dimension(:), allocatable u_g
Definition: fldvar_mod.f:87
param::top
integer top
Definition: param_mod.f:45
cutcell::area_w_cut
double precision, dimension(:), allocatable area_w_cut
Definition: cutcell_mod.f:134
cutcell::blocked_v_cell_at
logical, dimension(:), allocatable blocked_v_cell_at
Definition: cutcell_mod.f:365
cutcell::oneodz_t_v
double precision, dimension(:), allocatable oneodz_t_v
Definition: cutcell_mod.f:322
get_del_h
subroutine get_del_h(IJK, TYPE_OF_CELL, X0, Y0, Z0, Del_H, Nx, Ny, Nz)
Definition: get_delh.f:19
cutcell::z_w
double precision, dimension(:), allocatable z_w
Definition: cutcell_mod.f:61
geometry::vol
double precision, dimension(:), allocatable vol
Definition: geometry_mod.f:212
geometry::odz_t
double precision, dimension(:), allocatable odz_t
Definition: geometry_mod.f:125
bc::bc_ww_g
double precision, dimension(dimension_bc) bc_ww_g
Definition: bc_mod.f:319
param::bottom
integer bottom
Definition: param_mod.f:49
cutcell::y_u
double precision, dimension(:), allocatable y_u
Definition: cutcell_mod.f:50
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
bc
Definition: bc_mod.f:23
calc_tau_w_g
subroutine calc_tau_w_g(lTAU_W_G, lctau_w_g)
Definition: tau_w_g.f:62