MFIX  2016-1
tau_u_g.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: CALC_Tau_U_g C
4 ! Purpose: Cross terms in the gradient of stress in U_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 u-momentum equation that cannot be C
11 ! cast in the form: mu.grad(u). These components are stored in the C
12 ! passed variable, which is then used as a source of the u-momentum C
13 ! equation. C
14 ! C
15 ! The following u component viscous stress tensor terms are C
16 ! calculated here: C
17 ! > Combined part of 1/x d/dx (x.tau_xx) xdxdydz and C
18 ! -tau_zz/x xdxdydz => C
19 ! 1/x d/dx (x.lambda.trcD) xdxdydz - (lambda/x.trcD) xdxdydz => C
20 ! d/dx (lambda.trcD) xdxdydz => C
21 ! delta (lambda.trcD)Ap|E-W C
22 ! > part of 1/x d/dx(x.tau_xx) xdxdydz => C
23 ! 1/x d/dx (x.mu.du/dx) xdxdydz => C
24 ! delta (mu du/dx)Ayz |E-W C
25 ! > part of d/dy (tau_xy) xdxdydz => C
26 ! d/dy (mu.dv/dx) xdxdydz => C
27 ! delta (mu.dv/dx)Axz |N-S C
28 ! > part of 1/x d/dz (tau_xz) xdxdydz => C
29 ! 1/x d/dz (mu.dw/dx) xdxdydz => C
30 ! delta (mu.dw/dx)Axy |T-B C
31 ! CYLINDRICAL TERMS: C
32 ! > part of 1/x d/dz (tau_xz) xdxdydz => C
33 ! 1/x d/dz (mu.(-w/x)) xdxdydz => C
34 ! delta (mu(-w/x))Axy |T-B C
35 ! > part of -tau_zz/x xdxdydz => C
36 ! -(2.mu/x).(1/x).dw/dz xdxdydz C
37 ! delta (2.mu/x.1/x.dw/dz)Vp |E-W C
38 ! C
39 ! The following u-component CYLINDRICAL viscous stress tensor terms C
40 ! are not calculated here. They are handled in source_u_g: C
41 ! > part of the term -tau_zz/x xdxdydz => C
42 ! -2.mu/x. u/x xdxdydz C
43 ! delta(-2.mu.u/xp^2)Vp |p C
44 ! C
45 ! To reconstitute the complete u-momentum gas phase viscous stress C
46 ! tensor would require including the term calculated in source_u_g C
47 ! and the 'diffusional' components (i.e., those of the form C
48 ! mu.grad(u) C
49 ! C
50 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
51  SUBROUTINE calc_tau_u_g(lTAU_U_G, lctau_u_g)
52 
53 ! Modules
54 !---------------------------------------------------------------------//
55  USE param, only: dimension_3
56  USE param1, only: zero, half
57  USE constant
58  USE physprop
59  USE fldvar
60  USE visc_g
61  USE run
62  USE toleranc, only: dil_ep_s
63  USE geometry
64  USE indices
65  USE is
66  USE compar
67  USE sendrecv
68  USE functions
69  USE cutcell
70  IMPLICIT NONE
71 
72 ! Dummy arguments
73 !---------------------------------------------------------------------//
74 ! TAU_U_g
75  DOUBLE PRECISION, INTENT(OUT) :: lTAU_U_g(dimension_3)
76 ! cTAU_U_g
77  DOUBLE PRECISION, INTENT(OUT) :: lcTAU_U_g(dimension_3)
78 
79 ! Local variables
80 !---------------------------------------------------------------------//
81 ! Indices
82  INTEGER :: I, J, K, IP, JM, KM
83  INTEGER :: IJK, IJKE, IJKN, IJKS, IJKT, IJKB
84  INTEGER :: IJKNE, IJKSE, IJKTE, IJKBE
85  INTEGER :: IPJK, IMJK, IJMK, IJKM
86  INTEGER :: IPJMK, IPJKM
87 ! Average volume fraction
88  DOUBLE PRECISION :: EPGA
89 ! Average viscosity
90  DOUBLE PRECISION :: MU_gte, MU_gbe, MUGA
91 ! Average dW/Xdz
92  DOUBLE PRECISION :: dWoXdz
93 ! Source terms (Surface)
94  DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
95 ! Source terms (Volumetric)
96  DOUBLE PRECISION :: Vtzb
97 !---------------------------------------------------------------------//
98 
99  IF((.NOT.cartesian_grid).OR.(cg_safe_mode(3)==1)) THEN
100 
101 !$omp parallel do default(none) &
102 !$omp private(I, J, K, IP, JM, KM, IJK, &
103 !$omp IJKE, IJKN, IJKS, IJKT, IJKB, &
104 !$omp IJKNE, IJKSE, IJKTE, IJKBE, &
105 !$omp IPJK, IMJK, IJMK, IJKM, IPJMK, IPJKM, &
106 !$omp EPGA, MUGA, MU_GTE, MU_GBE, SBV, SSX, SSY, SSZ, &
107 !$omp DWOXDZ, VTZB) &
108 !$omp shared(ijkstart3, ijkend3, i_of, j_of, k_of, ip1, jm1, km1, &
109 !$omp do_k, cylindrical, ltau_u_g, lctau_u_g, &
110 !$omp ep_g, epmu_gt, lambda_gt, trd_g, u_g, v_g, w_g, &
111 !$omp axy_u, axz_u, ayz, ayz_u, vol_u, &
112 !$omp ox_e, odx_e, ox, odx, odz, eplambda_gt)
113  DO ijk = ijkstart3, ijkend3
114  i = i_of(ijk)
115  ijke = east_of(ijk)
116  epga = avg_x(ep_g(ijk),ep_g(ijke),i)
117  IF (.NOT.ip_at_e(ijk) .AND. epga>dil_ep_s) THEN
118  j = j_of(ijk)
119  k = k_of(ijk)
120  ip = ip1(i)
121  jm = jm1(j)
122  km = km1(k)
123 
124  ipjk = ip_of(ijk)
125  imjk = im_of(ijk)
126  ijmk = jm_of(ijk)
127  ijkm = km_of(ijk)
128  ipjmk = jm_of(ipjk)
129  ipjkm = ip_of(ijkm)
130 
131  ijkn = north_of(ijk)
132  ijkne = east_of(ijkn)
133  ijks = south_of(ijk)
134  ijkse = east_of(ijks)
135  ijkt = top_of(ijk)
136  ijkte = east_of(ijkt)
137  ijkb = bottom_of(ijk)
138  ijkbe = east_of(ijkb)
139 
140 
141 ! Surface forces at i+1/2, j, k
142 ! bulk viscosity term
143 ! combines part of 1/x d/dx (x.tau_xx) xdxdydz and -tau_zz/x xdxdydz =>
144 ! combines 1/x d/dx (x.lambda.trcD) xdxdydz - (lambda/x.trcD) xdxdydz =>
145 ! d/dx (lambda.trcD) xdxdydz
146 ! delta (lambda.trcD)Ap |E-W : at (i+1 - i-1), j, k
147  sbv = (eplambda_gt(ijke)*trd_g(ijke)-&
148  eplambda_gt(ijk)*trd_g(ijk))*ayz(ijk)
149 
150 ! shear stress terms at i+1/2, j, k
151 ! part of 1/x d/dx(x.tau_xx) xdxdydz =>
152 ! 1/x d/dx (x.mu.du/dx) xdxdydz =>
153 ! delta (mu du/dx)Ayz |E-W : at (i+1 - i-1), j, k
154  ssx = epmu_gt(ijke)*(u_g(ipjk)-u_g(ijk))*odx(ip)*&
155  ayz_u(ijk) - &
156  epmu_gt(ijk)*(u_g(ijk)-u_g(imjk))*odx(i)*&
157  ayz_u(imjk)
158 
159 ! part of d/dy (tau_xy) xdxdydz =>
160 ! d/dy (mu.dv/dx) xdxdydz =>
161 ! delta (mu.dv/dx)Axz |N-S : at i+1/2, (j+1/2 - j-1/2), k
162  ssy = avg_x_h(avg_y_h(epmu_gt(ijk),epmu_gt(ijkn),j),&
163  avg_y_h(epmu_gt(ijke),epmu_gt(ijkne),j),i)*&
164  (v_g(ipjk)-v_g(ijk))*odx_e(i)*axz_u(ijk) - &
165  avg_x_h(avg_y_h(epmu_gt(ijks),epmu_gt(ijk),jm),&
166  avg_y_h(epmu_gt(ijkse),epmu_gt(ijke),jm),i)*&
167  (v_g(ipjmk)-v_g(ijmk))*odx_e(i)*axz_u(ijmk)
168 
169 ! part of 1/x d/dz (tau_xz) xdxdydz =>
170 ! 1/x d/dz (mu.dw/dx) xdxdydz =>
171 ! delta (mu.dw/dx)Axy |T-B : at i+1/2, j, (k+1/2 - k-1/2)
172  mu_gte = avg_x_h(avg_z_h(epmu_gt(ijk),epmu_gt(ijkt),k),&
173  avg_z_h(epmu_gt(ijke),epmu_gt(ijkte),k),i)
174  mu_gbe = avg_x_h(avg_z_h(epmu_gt(ijkb),epmu_gt(ijk),km),&
175  avg_z_h(epmu_gt(ijkbe),epmu_gt(ijke),km),i)
176  ssz = mu_gte*(w_g(ipjk)-w_g(ijk))*odx_e(i)*&
177  axy_u(ijk) - &
178  mu_gbe*(w_g(ipjkm)-w_g(ijkm))*odx_e(i)*&
179  axy_u(ijkm)
180 
181 
182 ! Special terms for cylindrical coordinates
183  IF (cylindrical) THEN
184 ! part of 1/x d/dz (tau_xz) xdxdydz =>
185 ! 1/x d/dz (mu.(-w/x)) xdxdydz =>
186 ! delta (mu(-w/x))Axy |T-B : at i+1/2, j, (k+1/2- k-1/2)
187  ssz = ssz - ( &
188  mu_gte*(half*(w_g(ipjk)+w_g(ijk))*ox_e(i))*&
189  axy_u(ijk)-&
190  mu_gbe*(half*(w_g(ipjkm)+w_g(ijkm))*ox_e(i))*&
191  axy_u(ijkm) )
192 
193 ! part of -tau_zz/x xdxdydz =>
194 ! -(2.mu/x).(1/x).dw/dz xdxdydz
195 ! delta (2.mu/x.1/x.dw/dz)Vp |p : at i+1/2, j, k
196  muga = avg_x(epmu_gt(ijk),epmu_gt(ijke),i)
197  dwoxdz = half*((w_g(ijk)-w_g(ijkm))*ox(i)*odz(k)+&
198  (w_g(ipjk)-w_g(ipjkm))*ox(ip)*odz(k))
199  vtzb = -2.d0*muga*ox_e(i)*dwoxdz
200  ELSE
201  vtzb = zero
202  ENDIF
203 
204 ! Add the terms
205  ltau_u_g(ijk) = sbv + ssx + ssy + ssz + vtzb*vol_u(ijk)
206 
207 ! Also calculate and store the full gas phase viscous stress tensor
208  CALL get_full_tau_u_g(ijk, ltau_u_g, lctau_u_g)
209 
210  ELSE
211  ltau_u_g(ijk) = zero
212  lctau_u_g(ijk) = zero
213  ENDIF ! end if (.NOT.IP_AT_E(IJK) .AND. EPGA>DIL_EP_S)
214  ENDDO ! end do ijk
215 !$omp end parallel do
216 
217  ELSE
218 ! if cartesian grid
219  CALL calc_cg_tau_u_g(ltau_u_g, lctau_u_g)
220  ENDIF
221 
222  call send_recv(ltau_u_g,2)
223  call send_recv(lctau_u_g,2)
224 
225  RETURN
226 
227  CONTAINS
228 
229  include 'fun_avg.inc'
230 
231  END SUBROUTINE calc_tau_u_g
232 
233 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
234 ! C
235 ! Subroutine: CALC_CG_Tau_U_g C
236 ! Purpose: Cross terms in the gradient of stress in U_g momentum C
237 ! based on cartesian grid cut cell. C
238 ! C
239 ! Author: Jeff Dietiker Date: 01-Jul-09 C
240 ! C
241 ! C
242 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
243  SUBROUTINE calc_cg_tau_u_g(lTAu_U_G, lctau_u_g)
244 
245 ! Modules
246 !---------------------------------------------------------------------//
247 
248  USE bc
249  USE compar
250  USE constant
251  USE cutcell
252  USE fldvar
253  USE fun_avg
254  USE functions
255  USE geometry
256  USE indices
257  USE is
258  USE param, only: dimension_3
259  USE param1, only: zero, half, undefined
260  USE physprop
261  USE quadric
262  USE run
263  USE rxns
264  USE toleranc, only: dil_ep_s
265  USE visc_g
266 
267  IMPLICIT NONE
268 ! Dummy arguments
269 !---------------------------------------------------------------------//
270 ! TAU_U_g
271  DOUBLE PRECISION, INTENT(OUT) :: lTAU_U_g(dimension_3)
272 ! cTAU_U_g
273  DOUBLE PRECISION, INTENT(OUT) :: lcTAU_U_g(dimension_3)
274 
275 ! Local variables
276 !---------------------------------------------------------------------//
277 ! Indices
278  INTEGER :: I, J, K, IP, JM, KM
279  INTEGER :: IJK, IJKE, IJKN, IJKS, IJKT, IJKB
280  INTEGER :: IJKNE, IJKSE, IJKTE, IJKBE
281  INTEGER :: IPJK, IMJK, IJMK, IJKM
282  INTEGER :: IPJMK, IPJKM
283 ! Average volume fraction
284  DOUBLE PRECISION :: EPGA
285 ! Average viscosity
286  DOUBLE PRECISION :: MU_gte, MU_gbe
287 ! Source terms (Surface)
288  DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
289 ! Cartesian grid variables
290  DOUBLE PRECISION :: DEL_H, Nx, Ny, Nz
291  LOGICAL :: V_NODE_AT_NE, V_NODE_AT_NW, V_NODE_AT_SE, V_NODE_AT_SW
292  LOGICAL :: W_NODE_AT_TE, W_NODE_AT_TW, W_NODE_AT_BE, W_NODE_AT_BW
293  DOUBLE PRECISION :: dvdx_at_N, dvdx_at_S
294  DOUBLE PRECISION :: dwdx_at_T, dwdx_at_B
295  DOUBLE PRECISION :: Xi, Yi, Zi, Vi, Wi, Sx, Sy, Sz
296  DOUBLE PRECISION :: MU_GT_CUT, SSY_CUT, SSZ_CUT
297  DOUBLE PRECISION :: UW_g, VW_g, WW_g
298  INTEGER :: BCV
299  INTEGER :: BCT
300 !---------------------------------------------------------------------//
301 
302 !$omp parallel do default(none) &
303 !$omp private(I, J, K, JM, KM, IP, IJK, &
304 !$omp IJKE, IJKN, IJKS, IJKT, IJKB, &
305 !$omp IJKNE, IJKSE, IJKTE, IJKBE, &
306 !$omp IPJK, IMJK, IJMK, IJKM, IPJMK, IPJKM, &
307 !$omp EPGA, MU_GTE, MU_GBE, SBV, SSX, SSY, SSZ, &
308 !$omp BCV, BCT, BC_TYPE_ENUM, NOC_UG, uw_g, vw_g, ww_g, &
309 !$omp del_h, nx, ny, nz, xi, yi, zi, sx, sy, sz, vi, wi, &
310 !$omp v_node_at_ne, v_node_at_sw, v_node_at_se, v_node_at_nw, &
311 !$omp w_node_at_bw, w_node_at_be, w_node_at_tw, w_node_at_te, &
312 !$omp dvdx_at_S, dvdx_at_N, dwdx_at_B, dwdx_at_T, &
313 !$omp mu_gt_cut, cut_tau_ug, ssy_cut, ssz_cut) &
314 !$omp shared(ijkstart3, ijkend3, i_of, j_of, k_of, ip1, jm1, km1, &
315 !$omp do_k, ltau_u_g, lctau_u_g, &
316 !$omp ep_g, epmu_gt, lambda_gt, trd_g, u_g, v_g, w_g, &
317 !$omp axy_u, axz_u, ayz, ayz_u, vol, &
318 !$omp bc_type, bc_u_id, bc_uw_g, bc_vw_g, bc_ww_g, bc_hw_g, &
319 !$omp oneodx_e_u, oneodx_e_v, oneodx_e_w, &
320 !$omp x_u, y_u, z_u, x_v, y_v, z_v, x_w, y_w, z_w, &
321 !$omp wall_v_at, wall_w_at, area_u_cut, cut_u_cell_at, &
322 !$omp blocked_v_cell_at, blocked_w_cell_at, eplambda_gt)
323 
324  DO ijk = ijkstart3, ijkend3
325  i = i_of(ijk)
326  ijke = east_of(ijk)
327  epga = avg_x(ep_g(ijk),ep_g(ijke),i)
328  IF ( .NOT.ip_at_e(ijk) .AND. epga>dil_ep_s) THEN
329  j = j_of(ijk)
330  k = k_of(ijk)
331  ip = ip1(i)
332  jm = jm1(j)
333  km = km1(k)
334 
335  ipjk = ip_of(ijk)
336  imjk = im_of(ijk)
337  ijmk = jm_of(ijk)
338  ijkm = km_of(ijk)
339  ipjmk = jm_of(ipjk)
340  ipjkm = ip_of(ijkm)
341 
342  ijkn = north_of(ijk)
343  ijkne = east_of(ijkn)
344  ijks = south_of(ijk)
345  ijkse = east_of(ijks)
346  ijkt = top_of(ijk)
347  ijkte = east_of(ijkt)
348  ijkb = bottom_of(ijk)
349  ijkbe = east_of(ijkb)
350 
351 
352 ! bulk viscosity term
353  sbv = (eplambda_gt(ijke)*trd_g(ijke)) * ayz_u(ijk) - &
354  (eplambda_gt(ijk) *trd_g(ijk) ) * ayz_u(imjk)
355 
356 ! shear stress terms
357  IF(.NOT.cut_u_cell_at(ijk)) THEN
358  ssx = epmu_gt(ijke)*(u_g(ipjk)-u_g(ijk))*&
359  oneodx_e_u(ijk)*ayz_u(ijk) - &
360  epmu_gt(ijk)*(u_g(ijk)-u_g(imjk))*&
361  oneodx_e_u(imjk)*ayz_u(imjk)
362 
363  ssy = avg_x_h(avg_y_h(epmu_gt(ijk),epmu_gt(ijkn),j),&
364  avg_y_h(epmu_gt(ijke),epmu_gt(ijkne),j),i)*&
365  (v_g(ipjk)-v_g(ijk))*oneodx_e_v(ijk)*axz_u(ijk) - &
366  avg_x_h(avg_y_h(epmu_gt(ijks),epmu_gt(ijk),jm),&
367  avg_y_h(epmu_gt(ijkse),epmu_gt(ijke),jm),i)*&
368  (v_g(ipjmk)-v_g(ijmk))*oneodx_e_v(ijmk)*axz_u(ijmk)
369 
370  IF(do_k) THEN
371  mu_gte = avg_x_h(avg_z_h(epmu_gt(ijk),epmu_gt(ijkt),k),&
372  avg_z_h(epmu_gt(ijke),epmu_gt(ijkte),k),i)
373  mu_gbe = avg_x_h(avg_z_h(epmu_gt(ijkb),epmu_gt(ijk),km),&
374  avg_z_h(epmu_gt(ijkbe),epmu_gt(ijke),km),i)
375  ssz = mu_gte*(w_g(ipjk)-w_g(ijk))*&
376  oneodx_e_w(ijk)*axy_u(ijk) - &
377  mu_gbe*(w_g(ipjkm)-w_g(ijkm))*&
378  oneodx_e_w(ijkm)*axy_u(ijkm)
379  ELSE
380  ssz = zero
381  ENDIF
382 
383 ! cut cell modifications
384 !---------------------------------------------------------------------//
385  ELSE
386 
387  bcv = bc_u_id(ijk)
388  IF(bcv > 0 ) THEN
389  bct = bc_type_enum(bcv)
390  ELSE
391  bct = none
392  ENDIF
393 
394  SELECT CASE (bct)
395  CASE (cg_nsw,cg_mi)
396  cut_tau_ug = .true.
397  noc_ug = .true.
398  uw_g = zero
399  vw_g = zero
400  ww_g = zero
401  CASE (cg_fsw)
402  cut_tau_ug = .false.
403  noc_ug = .false.
404  uw_g = zero
405  vw_g = zero
406  ww_g = zero
407  CASE(cg_psw)
408  IF(bc_hw_g(bc_u_id(ijk))==undefined) THEN ! same as NSW
409  cut_tau_ug = .true.
410  noc_ug = .true.
411  uw_g = bc_uw_g(bcv)
412  vw_g = bc_vw_g(bcv)
413  ww_g = bc_ww_g(bcv)
414  ELSEIF(bc_hw_g(bc_u_id(ijk))==zero) THEN ! same as FSW
415  cut_tau_ug = .false.
416  noc_ug = .false.
417  ELSE ! partial slip
418  cut_tau_ug = .false.
419  noc_ug = .false.
420  uw_g = zero
421  vw_g = zero
422  ww_g = zero
423  ENDIF
424  CASE (none)
425  ltau_u_g(ijk) = zero
426  cycle
427  END SELECT
428 
429  IF(cut_tau_ug) THEN
430  mu_gt_cut = (vol(ijk)*epmu_gt(ijk) + &
431  vol(ipjk)*epmu_gt(ijke))/&
432  (vol(ijk) + vol(ipjk))
433  ELSE
434  mu_gt_cut = zero
435  ENDIF
436 
437 ! SSX:
438  CALL get_del_h(ijk,'U_MOMENTUM', x_u(ijk), &
439  y_u(ijk), z_u(ijk), del_h, nx, ny, nz)
440 
441  ssx = epmu_gt(ijke)*(u_g(ipjk)-u_g(ijk))*&
442  oneodx_e_u(ijk)*ayz_u(ijk) - &
443  epmu_gt(ijk)*(u_g(ijk)-u_g(imjk))*&
444  oneodx_e_u(imjk)*ayz_u(imjk) - &
445  mu_gt_cut * (u_g(ijk) - uw_g) / del_h * &
446  (nx**2) * area_u_cut(ijk)
447 
448 ! SSY:
449  v_node_at_ne = ((.NOT.blocked_v_cell_at(ipjk)).AND.&
450  (.NOT.wall_v_at(ipjk)))
451  v_node_at_nw = ((.NOT.blocked_v_cell_at(ijk)).AND.&
452  (.NOT.wall_v_at(ijk)))
453  v_node_at_se = ((.NOT.blocked_v_cell_at(ipjmk)).AND.&
454  (.NOT.wall_v_at(ipjmk)))
455  v_node_at_sw = ((.NOT.blocked_v_cell_at(ijmk)).AND.&
456  (.NOT.wall_v_at(ijmk)))
457 
458  IF(v_node_at_ne.AND.v_node_at_nw) THEN
459  vi = half * (v_g(ipjk) + v_g(ijk))
460  xi = half * (x_v(ipjk) + x_v(ijk))
461  yi = half * (y_v(ipjk) + y_v(ijk))
462  zi = half * (z_v(ipjk) + z_v(ijk))
463  sx = x_v(ipjk) - x_v(ijk)
464  sy = y_v(ipjk) - y_v(ijk)
465  sz = z_v(ipjk) - z_v(ijk)
466  CALL get_del_h(ijk,'U_MOMENTUM',xi, yi, zi, &
467  del_h, nx, ny, nz)
468 
469  dvdx_at_n = (v_g(ipjk)-v_g(ijk)) * oneodx_e_v(ijk)
470  IF(noc_ug) dvdx_at_n = dvdx_at_n - ( (vi-vw_g)*&
471  oneodx_e_v(ijk)/del_h * (sy*ny+sz*nz) )
472  ELSE
473  dvdx_at_n = zero
474  ENDIF
475 
476  IF(v_node_at_se.AND.v_node_at_sw) THEN
477  vi = half * (v_g(ipjmk) + v_g(ijmk))
478  xi = half * (x_v(ipjmk) + x_v(ijmk))
479  yi = half * (y_v(ipjmk) + y_v(ijmk))
480  zi = half * (z_v(ipjmk) + z_v(ijmk))
481  sx = x_v(ipjmk) - x_v(ijmk)
482  sy = y_v(ipjmk) - y_v(ijmk)
483  sz = z_v(ipjmk) - z_v(ijmk)
484  CALL get_del_h(ijk,'U_MOMENTUM', xi, yi, zi, &
485  del_h, nx, ny, nz)
486 
487  dvdx_at_s = (v_g(ipjmk)-v_g(ijmk)) * oneodx_e_v(ijmk)
488  IF(noc_ug) dvdx_at_s = dvdx_at_s - ( (vi-vw_g)*&
489  oneodx_e_v(ijmk)/del_h * (sy*ny+sz*nz) )
490  ELSE
491  dvdx_at_s = zero
492  ENDIF
493 
494  IF(v_node_at_nw) THEN
495  CALL get_del_h(ijk,'U_MOMENTUM', x_v(ijk), &
496  y_v(ijk), z_v(ijk), del_h, nx, ny, nz)
497  ssy_cut = -mu_gt_cut * (v_g(ijk)-vw_g)/del_h * &
498  (nx*ny) * area_u_cut(ijk)
499  ELSE
500  ssy_cut = zero
501  ENDIF
502 
503  ssy = avg_x_h(avg_y_h(epmu_gt(ijk),epmu_gt(ijkn),j),&
504  avg_y_h(epmu_gt(ijke),epmu_gt(ijkne),j),i)*&
505  dvdx_at_n*axz_u(ijk) - &
506  avg_x_h(avg_y_h(epmu_gt(ijks),epmu_gt(ijk),jm),&
507  avg_y_h(epmu_gt(ijkse),epmu_gt(ijke),jm),i)*&
508  dvdx_at_s*axz_u(ijmk) + ssy_cut
509 
510 ! SSZ:
511  IF(do_k) THEN
512  w_node_at_te = ((.NOT.blocked_w_cell_at(ipjk)).AND.&
513  (.NOT.wall_w_at(ipjk)))
514  w_node_at_tw = ((.NOT.blocked_w_cell_at(ijk)).AND.&
515  (.NOT.wall_w_at(ijk)))
516  w_node_at_be = ((.NOT.blocked_w_cell_at(ipjkm)).AND.&
517  (.NOT.wall_w_at(ipjkm)))
518  w_node_at_bw = ((.NOT.blocked_w_cell_at(ijkm)).AND.&
519  (.NOT.wall_w_at(ijkm)))
520 
521  IF(w_node_at_te.AND.w_node_at_tw) THEN
522  wi = half * (w_g(ipjk) + w_g(ijk))
523  xi = half * (x_w(ipjk) + x_w(ijk))
524  yi = half * (y_w(ipjk) + y_w(ijk))
525  zi = half * (z_w(ipjk) + z_w(ijk))
526  sx = x_w(ipjk) - x_w(ijk)
527  sy = y_w(ipjk) - y_w(ijk)
528  sz = z_w(ipjk) - z_w(ijk)
529  CALL get_del_h(ijk, 'U_MOMENTUM', xi, yi, zi, &
530  del_h, nx, ny, nz)
531 
532  dwdx_at_t = (w_g(ipjk)-w_g(ijk)) * oneodx_e_w(ijk)
533  IF(noc_ug) dwdx_at_t = dwdx_at_t - ( (wi-ww_g) * &
534  oneodx_e_w(ijk)/del_h * (sy*ny+sz*nz) )
535  ELSE
536  dwdx_at_t = zero
537  ENDIF
538 
539  IF(w_node_at_be.AND.w_node_at_bw) THEN
540  wi = half * (w_g(ipjkm) + w_g(ijkm))
541  xi = half * (x_w(ipjkm) + x_w(ijkm))
542  yi = half * (y_w(ipjkm) + y_w(ijkm))
543  zi = half * (z_w(ipjkm) + z_w(ijkm))
544  sx = x_w(ipjkm) - x_w(ijkm)
545  sy = y_w(ipjkm) - y_w(ijkm)
546  sz = z_w(ipjkm) - z_w(ijkm)
547 
548  CALL get_del_h(ijk,'U_MOMENTUM', xi, yi, zi,&
549  del_h, nx, ny, nz)
550 
551  dwdx_at_b = (w_g(ipjkm)-w_g(ijkm)) * oneodx_e_w(ijkm)
552  IF(noc_ug) dwdx_at_b = dwdx_at_b - ( (wi-ww_g) *&
553  oneodx_e_w(ijkm)/del_h * (sy*ny+sz*nz))
554  ELSE
555  dwdx_at_b = zero
556  ENDIF
557 
558  IF(w_node_at_tw) THEN
559  CALL get_del_h(ijk,'U_MOMENTUM', x_w(ijk), &
560  y_w(ijk), z_w(ijk), del_h, nx, ny, nz)
561  ssz_cut = -mu_gt_cut * (w_g(ijk)-ww_g)/del_h * &
562  (nx*nz) * area_u_cut(ijk)
563  ELSE
564  ssz_cut = zero
565  ENDIF
566 
567  mu_gte = avg_x_h(avg_z_h(epmu_gt(ijk),epmu_gt(ijkt),k),&
568  avg_z_h(epmu_gt(ijke),epmu_gt(ijkte),k),i)
569  mu_gbe = avg_x_h(avg_z_h(epmu_gt(ijkb),epmu_gt(ijk),km),&
570  avg_z_h(epmu_gt(ijkbe),epmu_gt(ijke),km),i)
571  ssz = mu_gte*dwdx_at_t*axy_u(ijk) - &
572  mu_gbe*dwdx_at_b*axy_u(ijkm) + ssz_cut
573  ELSE
574  ssz = zero
575  ENDIF ! DO_K
576 
577  ENDIF ! end if/else cut_u_cell_at
578 
579 ! Add the terms
580  ltau_u_g(ijk) = sbv + ssx + ssy + ssz
581 
582 ! Also calculate and store the full gas phase viscous stress tensor
583  CALL get_full_tau_u_g(ijk, ltau_u_g, lctau_u_g)
584  ELSE
585  ltau_u_g(ijk) = zero
586  lctau_u_g(ijk) = zero
587  ENDIF
588  ENDDO
589 !$omp end parallel do
590 
591  RETURN
592  END SUBROUTINE calc_cg_tau_u_g
593 
594 
595 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
596 ! C
597 ! Subroutine: GET_FULL_Tau_U_g C
598 ! Purpose: Calculate the divergence of the complete gas phase stress C
599 ! tensor including all terms in the u-direction. C
600 ! C
601 ! C
602 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
603  SUBROUTINE get_full_tau_u_g(IJK, ltau_u_g, lctau_u_g)
604 
605 ! Modules
606 !---------------------------------------------------------------------//
607  USE fldvar, only: u_g
608 
609  USE functions, only: east_of, flow_at_e
610  USE functions, only: im_of, jm_of, km_of
611  USE functions, only: ip_of, jp_of, kp_of
612 
613  USE fun_avg, only: avg_x
614 
615  USE geometry, only: ox_e, do_k, cylindrical, vol_u
616  USE indices, only: i_of
617 
618  USE param
619  USE param1, only: zero
620  USE visc_g, only: mu_gt, df_gu
621  IMPLICIT NONE
622 
623 ! Dummy arguments
624 !---------------------------------------------------------------------//
625 ! ijk index
626  INTEGER, INTENT(IN) :: IJK
627 ! TAU_U_g
628  DOUBLE PRECISION, INTENT(IN) :: lTAU_U_g(dimension_3)
629 ! cTAU_U_g
630  DOUBLE PRECISION, INTENT(INOUT) :: lcTAU_U_g(dimension_3)
631 
632 ! Local Variables
633 !---------------------------------------------------------------------//
634 ! indices
635  INTEGER :: I, IJKE
636  INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
637 ! average viscosity
638  DOUBLE PRECISION :: MUGA
639 ! source terms
640  DOUBLE PRECISION :: SSX, SSY, SSZ
641 ! cylindrical terms
642  DOUBLE PRECISION :: vtza
643 !---------------------------------------------------------------------//
644 
645  ipjk = ip_of(ijk)
646  ijpk = jp_of(ijk)
647  ijkp = kp_of(ijk)
648  imjk = im_of(ijk)
649  ijmk = jm_of(ijk)
650  ijkm = km_of(ijk)
651 
652 ! additional terms (found in source_u_g)
653  vtza = zero
654  IF (cylindrical) THEN
655  i = i_of(ijk)
656  ijke = east_of(ijk)
657 ! part of -tau_zz/x xdxdydz =>
658 ! -(2mu/x)*(u/x) xdxdydz =>
659 ! delta(-2.mu.u/x^2)V |p
660  muga = avg_x(mu_gt(ijk),mu_gt(ijke),i)
661  vtza = 2.d0*muga*ox_e(i)*ox_e(i)*vol_u(ijk)*u_g(ijk)
662  ENDIF
663 
664 ! convection terms: see conv_dif_u_g
665  ssx = zero
666  ssy = zero
667  ssz = zero
668  IF (flow_at_e(ijk)) THEN
669 ! part of 1/x d/dx (x.tau_xx) xdxdydz =>
670 ! 1/x d/dx (x.mu.du/dx) xdxdydz =>
671 ! delta (mu.du/dx.Ayz) |E-W : at (i+1 - i-1), j, k
672  ssx = df_gu(ijk,east)*(u_g(ipjk) - u_g(ijk)) - &
673  df_gu(ijk,west)*(u_g(ijk) - u_g(ijkm))
674 
675 ! part of d/dy (tau_xy) xdxdydz =>
676 ! d/dy (mu.du/dy) xdxdydz =>
677 ! delta (mu.du/dy.Axz) |N-S : at (i+1/2, j+1/2 - j-1/2, k)
678  ssy = df_gu(ijk,north)*(u_g(ijpk)-u_g(ijk)) - &
679  df_gu(ijk,south)*(u_g(ijk)-u_g(ijmk))
680 
681  IF (do_k) THEN
682 ! part of 1/x d/dz (tau_xz) xdxdydz =>
683 ! 1/x d/dz (mu/x.du/dz) xdxdydz =>
684 ! delta (mu/x.du/dz.Axy) |T-B : at (i+1/2, j, k+1/2 - k-1/2)
685  ssz = df_gu(ijk,top)*(u_g(ijkp)-u_g(ijk)) - &
686  df_gu(ijk,bottom)*(u_g(ijk)-u_g(ijkm))
687  ENDIF
688  ENDIF ! end if flow_at_e
689 
690 ! Add the terms
691  lctau_u_g(ijk) = (ltau_u_g(ijk) + ssx + ssy + ssz + vtza)
692 
693  RETURN
694  END SUBROUTINE get_full_tau_u_g
indices::ip1
integer, dimension(:), allocatable ip1
Definition: indices_mod.f:50
cutcell::y_v
double precision, dimension(:), allocatable y_v
Definition: cutcell_mod.f:55
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
bc::bc_uw_g
double precision, dimension(dimension_bc) bc_uw_g
Definition: bc_mod.f:313
cutcell::cut_u_cell_at
logical, dimension(:), allocatable cut_u_cell_at
Definition: cutcell_mod.f:356
cutcell::oneodx_e_w
double precision, dimension(:), allocatable oneodx_e_w
Definition: cutcell_mod.f:324
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
fldvar::ep_g
double precision, dimension(:), allocatable ep_g
Definition: fldvar_mod.f:15
geometry::ox_e
double precision, dimension(:), allocatable ox_e
Definition: geometry_mod.f:143
cutcell::oneodx_e_v
double precision, dimension(:), allocatable oneodx_e_v
Definition: cutcell_mod.f:320
geometry::odx
double precision, dimension(:), allocatable odx
Definition: geometry_mod.f:114
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
cutcell::x_u
double precision, dimension(:), allocatable x_u
Definition: cutcell_mod.f:49
rxns
Definition: rxns_mod.f:1
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
cutcell::noc_ug
logical noc_ug
Definition: cutcell_mod.f:389
bc::bc_type_enum
integer, dimension(dimension_bc) bc_type_enum
Definition: bc_mod.f:146
geometry::ayz_u
double precision, dimension(:), allocatable ayz_u
Definition: geometry_mod.f:218
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::ayz
double precision, dimension(:), allocatable ayz
Definition: geometry_mod.f:206
cutcell::z_v
double precision, dimension(:), allocatable z_v
Definition: cutcell_mod.f:56
is
Definition: is_mod.f:11
geometry::axz_u
double precision, dimension(:), allocatable axz_u
Definition: geometry_mod.f:220
cutcell::oneodx_e_u
double precision, dimension(:), allocatable oneodx_e_u
Definition: cutcell_mod.f:315
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
fun_avg
Definition: fun_avg_mod.f:1
cutcell::blocked_w_cell_at
logical, dimension(:), allocatable blocked_w_cell_at
Definition: cutcell_mod.f:366
sendrecv::send_recv
Definition: sendrecv_mod.f:75
cutcell::bc_u_id
integer, dimension(:), allocatable bc_u_id
Definition: cutcell_mod.f:434
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
cutcell::wall_w_at
logical, dimension(:), allocatable wall_w_at
Definition: cutcell_mod.f:128
visc_g::trd_g
double precision, dimension(:), allocatable trd_g
Definition: visc_g_mod.f:4
geometry::ox
double precision, dimension(:), allocatable ox
Definition: geometry_mod.f:140
cutcell::cut_tau_ug
logical cut_tau_ug
Definition: cutcell_mod.f:390
calc_tau_u_g
subroutine calc_tau_u_g(lTAU_U_G, lctau_u_g)
Definition: tau_u_g.f:52
visc_g::epmu_gt
double precision, dimension(:), allocatable epmu_gt
Definition: visc_g_mod.f:13
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
param::south
integer south
Definition: param_mod.f:41
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
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::do_k
logical do_k
Definition: geometry_mod.f:30
indices::km1
integer, dimension(:), allocatable km1
Definition: indices_mod.f:52
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
visc_g::df_gu
double precision, dimension(:,:), allocatable df_gu
Definition: visc_g_mod.f:31
fldvar::u_g
double precision, dimension(:), allocatable u_g
Definition: fldvar_mod.f:87
param::top
integer top
Definition: param_mod.f:45
geometry::vol_u
double precision, dimension(:), allocatable vol_u
Definition: geometry_mod.f:224
get_full_tau_u_g
subroutine get_full_tau_u_g(IJK, ltau_u_g, lctau_u_g)
Definition: tau_u_g.f:604
cutcell::area_u_cut
double precision, dimension(:), allocatable area_u_cut
Definition: cutcell_mod.f:132
cutcell::blocked_v_cell_at
logical, dimension(:), allocatable blocked_v_cell_at
Definition: cutcell_mod.f:365
calc_cg_tau_u_g
subroutine calc_cg_tau_u_g(lTAu_U_G, lctau_u_g)
Definition: tau_u_g.f:244
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::axy_u
double precision, dimension(:), allocatable axy_u
Definition: geometry_mod.f:222
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