MFIX  2016-1
tau_v_g.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: CALC_Tau_V_g C
4 ! Purpose: Cross terms in the gradient of stress in V_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 v-momentum equation that cannot be C
11 ! cast in the form: mu.grad(v). These components are stored in the C
12 ! passed variable, which is then used as a source of the v-momentum C
13 ! equation. C
14 ! C
15 ! The following v component viscous stress tensor terms are C
16 ! calculated here: C
17 ! > part of d/dy (tau_yy) xdxdydz => C
18 ! d/dy (lambda.trcD) xdxdydz => C
19 ! delta (lambda.trcD)Ap|N-S C
20 ! > part of 1/x d/dx(x.tau_xy) xdxdydz => C
21 ! 1/x d/dx (x.mu.du/dy) xdxdydz => C
22 ! delta (x.mu.du/dy)Ayz |E-W C
23 ! > part of d/dy (tau_xy) xdxdydz => C
24 ! d/dy (mu.dv/dy) xdxdydz => C
25 ! delta (mu.dv/dx)Axz |N-S C
26 ! > part of 1/x d/dz (tau_xz) xdxdydz => C
27 ! 1/x d/dz (mu.dw/dy) xdxdydz => C
28 ! delta (mu.dw/dx)Axy |T-B C
29 ! C
30 ! To reconstitute the full v-momentum gas phase viscous stress C
31 ! tensor would require including the the 'diffusional' components C
32 ! (i.e., those of the form mu.grad(v) C
33 ! C
34 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
35  SUBROUTINE calc_tau_v_g(lTAU_V_G, lctau_v_g)
36 
37 ! Modules
38 !---------------------------------------------------------------------//
39  USE param
40  USE param1
41  USE constant
42  USE physprop
43  USE fldvar
44  USE visc_g
45  USE run
46  USE toleranc
47  USE geometry
48  USE indices
49  USE is
50  USE sendrecv
51  USE compar
52  USE functions
53  USE cutcell
54  IMPLICIT NONE
55 
56 ! Dummy arguments
57 !---------------------------------------------------------------------//
58 ! TAU_V_g
59  DOUBLE PRECISION, INTENT(OUT) :: lTAU_V_g(dimension_3)
60 ! cTAU_V_g
61  DOUBLE PRECISION, INTENT(OUT) :: lcTAU_V_g(dimension_3)
62 
63 ! Local variables
64 !---------------------------------------------------------------------//
65 ! Indices
66  INTEGER :: I, J, K, IM, JP, KM
67  INTEGER :: IJK, IJKN, IJKE, IJKW, IJKT, IJKB
68  INTEGER :: IJKTN, IJKBN, IJKNE, IJKNW
69  INTEGER :: IJPK, IJMK, IMJK, IJKM
70  INTEGER :: IMJPK, IJPKM
71 ! Average volume fraction
72  DOUBLE PRECISION :: EPGA
73 ! Source terms (Surface)
74  DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
75 !---------------------------------------------------------------------//
76 
77  IF((.NOT.cartesian_grid).OR.(cg_safe_mode(4)==1)) THEN
78 
79 !$omp parallel do default(none) &
80 !$omp private(I, J, K, IJK, IM, JP, KM, &
81 !$omp IJKE, IJKW, IJKN, IJKT, IJKB, &
82 !$omp IJKTN, IJKBN, IJKNW, IJKNE, &
83 !$omp IJPK, IMJK, IJMK, IJKM, IJPKM, IMJPK, &
84 !$omp EPGA, SBV, SSX, SSY, SSZ) &
85 !$omp shared(ijkstart3, ijkend3, i_of, j_of, k_of, im1, jp1, km1, &
86 !$omp do_k, ltau_v_g, lctau_v_g, &
87 !$omp ep_g, lambda_gt, trd_g, epmu_gt, mu_g, u_g, v_g, w_g, &
88 !$omp ayz_v, axz_v, axz, axy_v, &
89 !$omp ody_n, ody)
90  DO ijk = ijkstart3, ijkend3
91  j = j_of(ijk)
92  ijkn = north_of(ijk)
93  epga = avg_y(ep_g(ijk),ep_g(ijkn),j)
94  IF ( .NOT.ip_at_n(ijk) .AND. epga>dil_ep_s) THEN
95  i = i_of(ijk)
96  k = k_of(ijk)
97  jp = jp1(j)
98  im = im1(i)
99  km = km1(k)
100 
101  ijpk = jp_of(ijk)
102  ijmk = jm_of(ijk)
103  imjk = im_of(ijk)
104  ijkm = km_of(ijk)
105  ijpkm = jp_of(ijkm)
106  imjpk = im_of(ijpk)
107 
108  ijke = east_of(ijk)
109  ijkne = east_of(ijkn)
110  ijkw = west_of(ijk)
111  ijknw = north_of(ijkw)
112  ijkt = top_of(ijk)
113  ijktn = north_of(ijkt)
114  ijkb = bottom_of(ijk)
115  ijkbn = north_of(ijkb)
116 
117 ! Surface forces at i, j+1/2, k
118 ! bulk viscosity term
119 ! part of d/dy (tau_yy) xdxdydz =>
120 ! d/dy (lambda.trcD) xdxdydz =>
121 ! delta (lambda.trcD)Ap|N-S : at (i, j+1 - j-1, k)
122  sbv = (lambda_gt(ijkn)*trd_g(ijkn)-&
123  lambda_gt(ijk)*trd_g(ijk))*axz(ijk)
124 
125 ! shear stress terms at i, j+1/2, k
126 
127 ! part of 1/x d/dx(x.tau_xy) xdxdydz =>
128 ! 1/x d/dx (x.mu.du/dy) xdxdydz =>
129 ! delta (x.mu.du/dy)Ayz |E-W : at (i+1/2 - i-1/2, j+1/2, k)
130  ssx = avg_y_h(avg_x_h(epmu_gt(ijk),epmu_gt(ijke),i),&
131  avg_x_h(epmu_gt(ijkn),epmu_gt(ijkne),i),j)*&
132  (u_g(ijpk)-u_g(ijk))*ody_n(j)*ayz_v(ijk) - &
133  avg_y_h(avg_x_h(epmu_gt(ijkw),epmu_gt(ijk),im),&
134  avg_x_h(epmu_gt(ijknw),epmu_gt(ijkn),im),j)*&
135  (u_g(imjpk)-u_g(imjk))*ody_n(j)*ayz_v(imjk)
136 
137 ! part of d/dy (tau_xy) xdxdydz =>
138 ! d/dy (mu.dv/dy) xdxdydz =>
139 ! delta (mu.dv/dx)Axz |N-S : at (i, j+1 - j-1, k)
140  ssy = epmu_gt(ijkn)*(v_g(ijpk)-v_g(ijk))*ody(jp)*&
141  axz_v(ijk) - &
142  epmu_gt(ijk)*(v_g(ijk)-v_g(ijmk))*ody(j)*&
143  axz_v(ijmk)
144 
145 ! part of 1/x d/dz (tau_xz) xdxdydz =>
146 ! 1/x d/dz (mu.dw/dy) xdxdydz =>
147 ! delta (mu.dw/dx)Axy |T-B : at (i, j+1/2, k+1/2 - k-1/2)
148  ssz = avg_y_h(avg_z_h(epmu_gt(ijk),epmu_gt(ijkt),k),&
149  avg_z_h(epmu_gt(ijkn),epmu_gt(ijktn),k),j)*&
150  (w_g(ijpk)-w_g(ijk))*ody_n(j)*axy_v(ijk) - &
151  avg_y_h(avg_z_h(epmu_gt(ijkb),epmu_gt(ijk),km),&
152  avg_z_h(epmu_gt(ijkbn),epmu_gt(ijkn),km),j)*&
153  (w_g(ijpkm)-w_g(ijkm))*ody_n(j)*axy_v(ijkm)
154 
155 ! Add the terms
156  ltau_v_g(ijk) = sbv + ssx + ssy + ssz
157 
158 ! Also calculate and store the full gas phase viscous stress tensor
159  CALL get_full_tau_v_g(ijk, ltau_v_g, lctau_v_g)
160 
161  ELSE
162  ltau_v_g(ijk) = zero
163  lctau_v_g(ijk) = zero
164  ENDIF ! end if (.NOT. IP_AT_N(IJK) .AND. EPGA>DIL_EP_S)
165  ENDDO ! end do ijk
166 !$omp end parallel do
167 
168  ELSE
169 ! cartesian grid case
170  CALL calc_cg_tau_v_g(ltau_v_g, lctau_v_g)
171  ENDIF
172 
173  call send_recv(ltau_v_g,2)
174  call send_recv(lctau_v_g,2)
175  RETURN
176 
177  CONTAINS
178 
179  include 'fun_avg.inc'
180 
181  END SUBROUTINE calc_tau_v_g
182 
183 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
184 ! C
185 ! Subroutine: CALC_CG_Tau_v_g C
186 ! Purpose: Cross terms in the gradient of stress in V_g momentum C
187 ! based on cartesian grid cut cell. C
188 ! C
189 ! Author: Jeff Dietiker Date: 01-Jul-09 C
190 ! C
191 ! C
192 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
193  SUBROUTINE calc_cg_tau_v_g(lTAU_V_G, lctau_v_g)
194 
195 ! Modules
196 !---------------------------------------------------------------------//
197  USE param
198  USE param1
199  USE constant
200  USE physprop
201  USE fldvar
202  USE visc_g
203  USE run
204  USE toleranc
205  USE geometry
206  USE indices
207  USE is
208  USE sendrecv
209  USE compar
210  USE fun_avg
211  USE functions
212 
213  USE cutcell
214 
215  USE bc
216  USE quadric
217  USE cutcell
218 
219  IMPLICIT NONE
220 
221 ! Dummy arguments
222 !---------------------------------------------------------------------//
223 ! TAU_V_g
224  DOUBLE PRECISION, INTENT(OUT) :: lTAU_V_g(dimension_3)
225 ! CTAU_V_g
226  DOUBLE PRECISION, INTENT(OUT) :: lcTAU_V_g(dimension_3)
227 
228 ! Local variables
229 !---------------------------------------------------------------------//
230 ! Indices
231  INTEGER :: I, J, K, IM, JP, KM
232  INTEGER :: IJK, IJKN, IJKE, IJKW, IJKT, IJKB
233  INTEGER :: IJKTN, IJKBN, IJKNE, IJKNW
234  INTEGER :: IJPK, IJMK, IMJK, IJKM
235  INTEGER :: IMJPK, IJPKM
236 ! Average volume fraction
237  DOUBLE PRECISION EPGA
238 ! Source terms (Surface)
239  DOUBLE PRECISION Sbv, Ssx, Ssy, Ssz
240 
241 ! Cartesian grid variables
242  DOUBLE PRECISION :: DEL_H, Nx, Ny, Nz
243  LOGICAL :: U_NODE_AT_NE, U_NODE_AT_NW, U_NODE_AT_SE, U_NODE_AT_SW
244  LOGICAL :: W_NODE_AT_TN, W_NODE_AT_TS, W_NODE_AT_BN, W_NODE_AT_BS
245  DOUBLE PRECISION :: dudy_at_E, dudy_at_W
246  DOUBLE PRECISION :: dwdy_at_T, dwdy_at_B
247  DOUBLE PRECISION :: Xi, Yi, Zi, Ui, Wi, Sx, Sy, Sz
248  DOUBLE PRECISION :: MU_GT_CUT, SSX_CUT, SSZ_CUT
249  DOUBLE PRECISION :: UW_g, VW_g, WW_g
250  INTEGER :: BCV
251  INTEGER :: BCT
252 !---------------------------------------------------------------------//
253 
254 !$omp parallel do default(none) &
255 !$omp private(I, J, K, IM, JP, KM, IJK, &
256 !$omp IJKE, IJKW, IJKN, IJKT, IJKB, &
257 !$omp IJKTN, IJKBN, IJKNW, IJKNE, &
258 !$omp IJPK, IMJK, IJMK, IJKM, IJPKM, IMJPK, &
259 !$omp EPGA, SBV, SSX, SSY, SSZ, &
260 !$omp BCV, BCT, BC_TYPE_ENUM, NOC_VG, uw_g, vw_g, ww_g, &
261 !$omp del_h, nx, ny, nz, xi, yi, zi, sx, sy, sz, wi, ui, &
262 !$omp u_node_at_sw, u_node_at_se, u_node_at_nw, u_node_at_ne, &
263 !$omp w_node_at_bs, w_node_at_bn, w_node_at_ts, w_node_at_tn, &
264 !$omp dudy_at_W, dudy_at_E, dwdy_at_B, dwdy_at_T, &
265 !$omp cut_tau_vg, mu_gt_cut, ssz_cut, ssx_cut) &
266 !$omp shared(ijkstart3, ijkend3, i_of, j_of, k_of, do_k, &
267 !$omp im1, jm1, jp1, km1, ltau_v_g, lctau_v_g, &
268 !$omp epmu_gt, ep_g, lambda_gt, trd_g, v_g, w_g, u_g, &
269 !$omp ayz_v, axz_v, axz, axy_v, vol, &
270 !$omp bc_type, bc_v_id, bc_hw_g, bc_uw_g, bc_vw_g, bc_ww_g, &
271 !$omp x_u, y_u, z_u, x_v, y_v, z_v, x_w, y_w, z_w, &
272 !$omp oneody_n_v, oneody_n_w, oneody_n_u, &
273 !$omp cut_v_cell_at, wall_u_at, wall_w_at, area_v_cut, &
274 !$omp blocked_u_cell_at, blocked_w_cell_at)
275 
276  DO ijk = ijkstart3, ijkend3
277  j = j_of(ijk)
278  ijkn = north_of(ijk)
279  epga = avg_y(ep_g(ijk),ep_g(ijkn),j)
280  IF ( .NOT.ip_at_n(ijk) .AND. epga>dil_ep_s) THEN
281  i = i_of(ijk)
282  k = k_of(ijk)
283  jp = jp1(j)
284  im = im1(i)
285  km = km1(k)
286 
287  ijpk = jp_of(ijk)
288  ijmk = jm_of(ijk)
289  imjk = im_of(ijk)
290  ijkm = km_of(ijk)
291  imjpk = im_of(ijpk)
292  ijpkm = jp_of(ijkm)
293 
294  ijke = east_of(ijk)
295  ijkne = east_of(ijkn)
296  ijkw = west_of(ijk)
297  ijknw = north_of(ijkw)
298  ijkt = top_of(ijk)
299  ijktn = north_of(ijkt)
300  ijkb = bottom_of(ijk)
301  ijkbn = north_of(ijkb)
302 
303 
304 ! bulk viscosity term
305  sbv = (lambda_gt(ijkn)*trd_g(ijkn)) * axz_v(ijk) - &
306  (lambda_gt(ijk) *trd_g(ijk) ) * axz_v(ijmk)
307 
308 ! shear stress terms
309  IF(.NOT.cut_v_cell_at(ijk)) THEN
310  ssx = avg_y_h(avg_x_h(epmu_gt(ijk),epmu_gt(ijke),i),&
311  avg_x_h(epmu_gt(ijkn),epmu_gt(ijkne),i),j)*&
312  (u_g(ijpk)-u_g(ijk))*oneody_n_u(ijk)*ayz_v(ijk) - &
313  avg_y_h(avg_x_h(epmu_gt(ijkw),epmu_gt(ijk),im),&
314  avg_x_h(epmu_gt(ijknw),epmu_gt(ijkn),im),j)*&
315  (u_g(imjpk)-u_g(imjk))*oneody_n_u(imjk)*ayz_v(imjk)
316 
317  ssy = epmu_gt(ijkn)*(v_g(ijpk)-v_g(ijk))*&
318  oneody_n_v(ijk)*axz_v(ijk) - &
319  epmu_gt(ijk)*(v_g(ijk)-v_g(ijmk))*&
320  oneody_n_v(ijmk)*axz_v(ijmk)
321 
322  IF(do_k) THEN
323  ssz = avg_y_h(avg_z_h(epmu_gt(ijk),epmu_gt(ijkt),k),&
324  avg_z_h(epmu_gt(ijkn),epmu_gt(ijktn),k),j)*&
325  (w_g(ijpk)-w_g(ijk))*oneody_n_w(ijk)*axy_v(ijk) - &
326  avg_y_h(avg_z_h(epmu_gt(ijkb),epmu_gt(ijk),km),&
327  avg_z_h(epmu_gt(ijkbn),epmu_gt(ijkn),km),j)*&
328  (w_g(ijpkm)-w_g(ijkm))*oneody_n_w(ijkm)*axy_v(ijkm)
329  ELSE
330  ssz = zero
331  ENDIF
332 
333 ! cut cell modifications
334 !---------------------------------------------------------------------//
335  ELSE
336 
337  bcv = bc_v_id(ijk)
338  IF(bcv > 0 ) THEN
339  bct = bc_type_enum(bcv)
340  ELSE
341  bct = none
342  ENDIF
343 
344  SELECT CASE (bct)
345  CASE (cg_nsw,cg_mi)
346  cut_tau_vg = .true.
347  noc_vg = .true.
348  uw_g = zero
349  vw_g = zero
350  ww_g = zero
351  CASE (cg_fsw)
352  cut_tau_vg = .false.
353  noc_vg = .false.
354  uw_g = zero
355  vw_g = zero
356  ww_g = zero
357  CASE(cg_psw)
358  IF(bc_hw_g(bc_v_id(ijk))==undefined) THEN ! same as NSW
359  cut_tau_vg = .true.
360  noc_vg = .true.
361  uw_g = bc_uw_g(bcv)
362  vw_g = bc_vw_g(bcv)
363  ww_g = bc_ww_g(bcv)
364  ELSEIF(bc_hw_g(bc_v_id(ijk))==zero) THEN ! same as FSW
365  cut_tau_vg = .false.
366  noc_vg = .false.
367  uw_g = zero
368  vw_g = zero
369  ww_g = zero
370  ELSE ! partial slip
371  cut_tau_vg = .false.
372  noc_vg = .false.
373  ENDIF
374  CASE (none)
375  ltau_v_g(ijk) = zero
376  cycle
377  END SELECT
378 
379 
380  IF(cut_tau_vg) THEN
381  mu_gt_cut = (vol(ijk)*epmu_gt(ijk) + &
382  vol(ijpk)*epmu_gt(ijkn))/&
383  (vol(ijk) + vol(ijpk))
384  ELSE
385  mu_gt_cut = zero
386  ENDIF
387 
388 ! SSX:
389  u_node_at_ne = ((.NOT.blocked_u_cell_at(ijpk)).AND.&
390  (.NOT.wall_u_at(ijpk)))
391  u_node_at_se = ((.NOT.blocked_u_cell_at(ijk)).AND.&
392  (.NOT.wall_u_at(ijk)))
393  u_node_at_nw = ((.NOT.blocked_u_cell_at(imjpk)).AND.&
394  (.NOT.wall_u_at(imjpk)))
395  u_node_at_sw = ((.NOT.blocked_u_cell_at(imjk)).AND.&
396  (.NOT.wall_u_at(imjk)))
397 
398  IF(u_node_at_ne.AND.u_node_at_se) THEN
399  ui = half * (u_g(ijpk) + u_g(ijk))
400  xi = half * (x_u(ijpk) + x_u(ijk))
401  yi = half * (y_u(ijpk) + y_u(ijk))
402  zi = half * (z_u(ijpk) + z_u(ijk))
403  sx = x_u(ijpk) - x_u(ijk)
404  sy = y_u(ijpk) - y_u(ijk)
405  sz = z_u(ijpk) - z_u(ijk)
406  CALL get_del_h(ijk, 'V_MOMENTUM', xi, yi, zi, &
407  del_h, nx, ny, nz)
408 
409  dudy_at_e = (u_g(ijpk) - u_g(ijk)) * oneody_n_u(ijk)
410  IF(noc_vg) dudy_at_e = dudy_at_e - ((ui-uw_g) * &
411  oneody_n_u(ijk)/del_h * (sx*nx+sz*nz) )
412  ELSE
413  dudy_at_e = zero
414  ENDIF
415 
416  IF(u_node_at_nw.AND.u_node_at_sw) THEN
417  ui = half * (u_g(imjpk) + u_g(imjk))
418  xi = half * (x_u(imjpk) + x_u(imjk))
419  yi = half * (y_u(imjpk) + y_u(imjk))
420  zi = half * (z_u(imjpk) + z_u(imjk))
421  sx = x_u(imjpk) - x_u(imjk)
422  sy = y_u(imjpk) - y_u(imjk)
423  sz = z_u(imjpk) - z_u(imjk)
424  CALL get_del_h(ijk, 'V_MOMENTUM', xi, yi, zi, &
425  del_h, nx, ny, nz)
426 
427  dudy_at_w = (u_g(imjpk) - u_g(imjk)) * oneody_n_u(imjk)
428  IF(noc_vg) dudy_at_w = dudy_at_w - ((ui-uw_g) * &
429  oneody_n_u(imjk)/del_h * (sx*nx+sz*nz) )
430  ELSE
431  dudy_at_w = zero
432  ENDIF
433 
434  IF(u_node_at_se) THEN
435  CALL get_del_h(ijk,'V_MOMENTUM', x_u(ijk), y_u(ijk), &
436  z_u(ijk), del_h, nx, ny, nz)
437  ssx_cut = - mu_gt_cut * (u_g(ijk) - uw_g) / del_h * &
438  (ny*nx) * area_v_cut(ijk)
439  ELSE
440  ssx_cut = zero
441  ENDIF
442 
443  ssx = avg_y_h(avg_x_h(epmu_gt(ijk),epmu_gt(ijke),i),&
444  avg_x_h(epmu_gt(ijkn),epmu_gt(ijkne),i),j)*&
445  dudy_at_e*ayz_v(ijk) - &
446  avg_y_h(avg_x_h(epmu_gt(ijkw),epmu_gt(ijk),im),&
447  avg_x_h(epmu_gt(ijknw),epmu_gt(ijkn),im),j)*&
448  dudy_at_w*ayz_v(imjk) + ssx_cut
449 
450 ! SSY:
451  CALL get_del_h(ijk, 'V_MOMENTUM', x_v(ijk), y_v(ijk), &
452  z_v(ijk), del_h, nx, ny, nz)
453 
454  ssy = epmu_gt(ijkn)*(v_g(ijpk)-v_g(ijk))*&
455  oneody_n_v(ijk)*axz_v(ijk) - &
456  epmu_gt(ijk)*(v_g(ijk)-v_g(ijmk))*&
457  oneody_n_v(ijmk)*axz_v(ijmk) - &
458  mu_gt_cut * (v_g(ijk) - vw_g) / del_h * &
459  (ny**2) * area_v_cut(ijk)
460 
461 ! SSZ:
462  IF(do_k) THEN
463  w_node_at_tn = ((.NOT.blocked_w_cell_at(ijpk)).AND.&
464  (.NOT.wall_w_at(ijpk)))
465  w_node_at_ts = ((.NOT.blocked_w_cell_at(ijk)).AND.&
466  (.NOT.wall_w_at(ijk)))
467  w_node_at_bn = ((.NOT.blocked_w_cell_at(ijpkm)).AND.&
468  (.NOT.wall_w_at(ijpkm)))
469  w_node_at_bs = ((.NOT.blocked_w_cell_at(ijkm)).AND.&
470  (.NOT.wall_w_at(ijkm)))
471 
472  IF(w_node_at_tn.AND.w_node_at_ts) THEN
473  wi = half * (w_g(ijpk) + w_g(ijk))
474  xi = half * (x_w(ijpk) + x_w(ijk))
475  yi = half * (y_w(ijpk) + y_w(ijk))
476  zi = half * (z_w(ijpk) + z_w(ijk))
477  sx = x_w(ijpk) - x_w(ijk)
478  sy = y_w(ijpk) - y_w(ijk)
479  sz = z_w(ijpk) - z_w(ijk)
480  CALL get_del_h(ijk, 'V_MOMENTUM', xi, yi, zi, &
481  del_h, nx, ny, nz)
482 
483  dwdy_at_t = (w_g(ijpk)-w_g(ijk)) * &
484  oneody_n_w(ijk)
485  IF(noc_vg) dwdy_at_t = dwdy_at_t - ((wi-ww_g) * &
486  oneody_n_w(ijk)/del_h*(sx*nx+sz*nz))
487  ELSE
488  dwdy_at_t = zero
489  ENDIF
490 
491  IF(w_node_at_bn.AND.w_node_at_bs) THEN
492  wi = half * (w_g(ijpkm) + w_g(ijkm))
493  xi = half * (x_w(ijpkm) + x_w(ijkm))
494  yi = half * (y_w(ijpkm) + y_w(ijkm))
495  zi = half * (z_w(ijpkm) + z_w(ijkm))
496  sx = x_w(ijpkm) - x_w(ijkm)
497  sy = y_w(ijpkm) - y_w(ijkm)
498  sz = z_w(ijpkm) - z_w(ijkm)
499  CALL get_del_h(ijk, 'V_MOMENTUM', xi, yi, zi,&
500  del_h, nx, ny, nz)
501 
502  dwdy_at_b = (w_g(ijpkm) - w_g(ijkm)) * &
503  oneody_n_w(ijkm)
504  IF(noc_vg) dwdy_at_b = dwdy_at_b - ((wi-ww_g) * &
505  oneody_n_w(ijkm)/del_h*(sx*nx+sz*nz))
506  ELSE
507  dwdy_at_b = zero
508  ENDIF
509 
510  IF(w_node_at_ts) THEN
511  CALL get_del_h(ijk, 'V_MOMENTUM', x_w(ijk), &
512  y_w(ijk), z_w(ijk), del_h, nx, ny, nz)
513  ssz_cut = - mu_gt_cut * (w_g(ijk) - ww_g) / &
514  del_h * (ny*nz) * area_v_cut(ijk)
515  ELSE
516  ssz_cut = zero
517  ENDIF
518 
519  ssz = avg_y_h(avg_z_h(epmu_gt(ijk),epmu_gt(ijkt),k),&
520  avg_z_h(epmu_gt(ijkn),epmu_gt(ijktn),k),j)*&
521  dwdy_at_t*axy_v(ijk) - &
522  avg_y_h(avg_z_h(epmu_gt(ijkb),epmu_gt(ijk),km),&
523  avg_z_h(epmu_gt(ijkbn),epmu_gt(ijkn),km),j)*&
524  dwdy_at_b*axy_v(ijkm) + ssz_cut
525  ELSE
526  ssz = zero
527  ENDIF ! DO_K
528 
529  ENDIF ! end if/else cut_v_cell_at
530 
531 ! Add the terms
532  ltau_v_g(ijk) = sbv + ssx + ssy + ssz
533 
534 ! Also calculate and store the full gas phase viscous stress tensor
535  CALL get_full_tau_v_g(ijk, ltau_v_g, lctau_v_g)
536  ELSE
537  ltau_v_g(ijk) = zero
538  lctau_v_g(ijk) = zero
539  ENDIF ! end if/else (.NOT. IP_AT_N(IJK) .AND. EPGA>DIL_EP_S)
540  ENDDO ! end do ijk
541 !$omp end parallel do
542 
543  RETURN
544  END SUBROUTINE calc_cg_tau_v_g
545 
546 
547 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
548 ! C
549 ! Subroutine: GET_FULL_Tau_V_g C
550 ! Purpose: Calculate the divergence of complete gas phase stress C
551 ! tensor including all terms in v-direction C
552 ! C
553 ! C
554 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
555  SUBROUTINE get_full_tau_v_g(IJK, ltau_v_g, lctau_v_g)
556 
557 ! Modules
558 !---------------------------------------------------------------------//
559  USE fldvar, only: v_g
560 
561  USE functions, only: flow_at_n
562  USE functions, only: im_of, jm_of, km_of
563  USE functions, only: ip_of, jp_of, kp_of
564 
565  USE geometry, only: do_k
566 
567  USE param
568  USE param1, only: zero
569  USE visc_g, only: df_gv
570  IMPLICIT NONE
571 
572 ! Dummy arguments
573 !---------------------------------------------------------------------//
574 ! ijk index
575  INTEGER, INTENT(IN) :: IJK
576 ! TAU_V_g
577  DOUBLE PRECISION, INTENT(IN) :: lTAU_V_g(dimension_3)
578 ! cTAU_V_g
579  DOUBLE PRECISION, INTENT(INOUT) :: lcTAU_V_g(dimension_3)
580 
581 ! Local Variables
582 !---------------------------------------------------------------------//
583 ! indices
584  INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
585 ! source terms
586  DOUBLE PRECISION :: SSX, SSY, SSZ
587 !---------------------------------------------------------------------//
588 
589  ipjk = ip_of(ijk)
590  ijpk = jp_of(ijk)
591  ijkp = kp_of(ijk)
592  imjk = im_of(ijk)
593  ijmk = jm_of(ijk)
594  ijkm = km_of(ijk)
595 
596 ! convection terms: see conv_dif_v_g
597  ssx = zero
598  ssy = zero
599  ssz = zero
600  IF (flow_at_n(ijk)) THEN
601 ! part of 1/x d/dx (x.tau_xx) xdxdydz =>
602 ! 1/x d/dx (x.mu.dv/dx) xdxdydz =>
603 ! delta (mu.dv/dx.Ayz) |E-W : at (i+1/2 - i-1/2), j+1/2, k
604  ssx = df_gv(ijk,east)*(v_g(ipjk) - v_g(ijk)) - &
605  df_gv(ijk,west)*(v_g(ijk) - v_g(ijkm))
606 
607 ! part of d/dy (tau_xy) xdxdydz =>
608 ! d/dy (mu.dv/dy) xdxdydz =>
609 ! delta (mu.dv/dy.Axz) |N-S : at (i, j+1 - j-1, k)
610  ssy = df_gv(ijk,north)*(v_g(ijpk)-v_g(ijk)) - &
611  df_gv(ijk,south)*(v_g(ijk)-v_g(ijmk))
612 
613  IF (do_k) THEN
614 ! part of 1/x d/dz (tau_xz) xdxdydz =>
615 ! 1/x d/dz (mu/x.dv/dz) xdxdydz =>
616 ! delta (mu/x.dv/dz.Axy) |T-B : at (i, j+1/2, k+1/2 - k-1/2)
617  ssz = df_gv(ijk,top)*(v_g(ijkp)-v_g(ijk)) - &
618  df_gv(ijk,bottom)*(v_g(ijk)-v_g(ijkm))
619  ENDIF
620  ENDIF ! end if flow_at_n
621 
622 ! Add the terms
623  lctau_v_g(ijk) = (ltau_v_g(ijk) + ssx + ssy + ssz)
624 
625  RETURN
626  END SUBROUTINE get_full_tau_v_g
627 
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
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
cutcell::noc_vg
logical noc_vg
Definition: cutcell_mod.f:389
fldvar::ep_g
double precision, dimension(:), allocatable ep_g
Definition: fldvar_mod.f:15
geometry::ody
double precision, dimension(:), allocatable ody
Definition: geometry_mod.f:116
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
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
visc_g::df_gv
double precision, dimension(:,:), allocatable df_gv
Definition: visc_g_mod.f:32
cutcell::cg_safe_mode
integer, dimension(10) cg_safe_mode
Definition: cutcell_mod.f:415
bc::bc_type_enum
integer, dimension(dimension_bc) bc_type_enum
Definition: bc_mod.f:146
get_full_tau_v_g
subroutine get_full_tau_v_g(IJK, ltau_v_g, lctau_v_g)
Definition: tau_v_g.f:556
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
cutcell::z_v
double precision, dimension(:), allocatable z_v
Definition: cutcell_mod.f:56
cutcell::oneody_n_u
double precision, dimension(:), allocatable oneody_n_u
Definition: cutcell_mod.f:316
is
Definition: is_mod.f:11
geometry::ayz_v
double precision, dimension(:), allocatable ayz_v
Definition: geometry_mod.f:227
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
geometry::ody_n
double precision, dimension(:), allocatable ody_n
Definition: geometry_mod.f:123
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::x_w
double precision, dimension(:), allocatable x_w
Definition: cutcell_mod.f:59
cutcell::area_v_cut
double precision, dimension(:), allocatable area_v_cut
Definition: cutcell_mod.f:133
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
cutcell::wall_w_at
logical, dimension(:), allocatable wall_w_at
Definition: cutcell_mod.f:128
geometry::axy_v
double precision, dimension(:), allocatable axy_v
Definition: geometry_mod.f:231
visc_g::trd_g
double precision, dimension(:), allocatable trd_g
Definition: visc_g_mod.f:4
indices::jp1
integer, dimension(:), allocatable jp1
Definition: indices_mod.f:51
cutcell::oneody_n_v
double precision, dimension(:), allocatable oneody_n_v
Definition: cutcell_mod.f:321
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
calc_cg_tau_v_g
subroutine calc_cg_tau_v_g(lTAU_V_G, lctau_v_g)
Definition: tau_v_g.f:194
fldvar::v_g
double precision, dimension(:), allocatable v_g
Definition: fldvar_mod.f:99
param::south
integer south
Definition: param_mod.f:41
cutcell::cut_tau_vg
logical cut_tau_vg
Definition: cutcell_mod.f:390
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
calc_tau_v_g
subroutine calc_tau_v_g(lTAU_V_G, lctau_v_g)
Definition: tau_v_g.f:36
geometry::axz
double precision, dimension(:), allocatable axz
Definition: geometry_mod.f:208
cutcell::bc_v_id
integer, dimension(:), allocatable bc_v_id
Definition: cutcell_mod.f:435
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
cutcell::cut_v_cell_at
logical, dimension(:), allocatable cut_v_cell_at
Definition: cutcell_mod.f:357
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
indices::km1
integer, dimension(:), allocatable km1
Definition: indices_mod.f:52
bc::bc_vw_g
double precision, dimension(dimension_bc) bc_vw_g
Definition: bc_mod.f:316
visc_g::lambda_gt
double precision, dimension(:), allocatable lambda_gt
Definition: visc_g_mod.f:17
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
fldvar::u_g
double precision, dimension(:), allocatable u_g
Definition: fldvar_mod.f:87
param::top
integer top
Definition: param_mod.f:45
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
cutcell::oneody_n_w
double precision, dimension(:), allocatable oneody_n_w
Definition: cutcell_mod.f:325
geometry::vol
double precision, dimension(:), allocatable vol
Definition: geometry_mod.f:212
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
geometry::axz_v
double precision, dimension(:), allocatable axz_v
Definition: geometry_mod.f:229
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