MFIX  2016-1
calc_trd_g.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: CALC_trD_g C
4 ! Purpose: Calculate the trace of the gas phase rate of strain C
5 ! tensor at i, j, k C
6 ! C
7 ! Author: M. Syamlal Date: 19-DEC-96 C
8 ! C
9 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
10  SUBROUTINE calc_trd_g(lTRD_G)
11 
12 ! Modules
13 !-----------------------------------------------
14  USE param, only: dimension_3
15  USE param1, only: zero
16  USE parallel
17  USE geometry
18  USE fldvar
19  USE indices
20  USE compar
21  USE sendrecv
22  USE functions
23  USE cutcell
24  IMPLICIT NONE
25 
26 ! Dummy arguments
27 !-----------------------------------------------
28 ! Strain rate tensor components for mth solids phase
29  DOUBLE PRECISION, INTENT(OUT) :: ltrD_g(dimension_3)
30 
31 ! Local variables
32 !-----------------------------------------------
33 ! Indices
34  INTEGER :: I, J, K, IJK, IMJK, IJMK, IJKM, IM
35 
36 !-----------------------------------------------
37 
38 
39 !!$omp parallel do private( IJK, I, J, K, IM, IMJK, IJMK, IJKM ) &
40 !!$omp& schedule(dynamic,chunk_size)
41  DO ijk = ijkstart3, ijkend3
42  IF (.NOT.wall_at(ijk)) THEN
43  i = i_of(ijk)
44  j = j_of(ijk)
45  k = k_of(ijk)
46  im = im1(i)
47  imjk = im_of(ijk)
48  ijmk = jm_of(ijk)
49  ijkm = km_of(ijk)
50 
51  IF (.NOT. cut_cell_at(ijk)) THEN
52 ! at i, j, k:
53 ! du/dx + dv/dy + 1/x dw/dz + u/x =
54 ! 1/x d(xu)/dx + dv/dy + 1/x dw/dz =
55  ltrd_g(ijk) = &
56  (x_e(i)*u_g(ijk)-x_e(im)*u_g(imjk))*ox(i)*odx(i) +&
57  (v_g(ijk)-v_g(ijmk))*ody(j) + &
58  (w_g(ijk)-w_g(ijkm))*(ox(i)*odz(k))
59  ELSE
60  CALL calc_cg_trd_g(ijk,ltrd_g)
61  ENDIF
62 
63  ELSE
64  ltrd_g(ijk) = zero
65  ENDIF
66  ENDDO
67 
68  RETURN
69  END SUBROUTINE calc_trd_g
70 
71 
72 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
73 ! C
74 ! Subroutine: CALC_CG_trD_g C
75 ! Purpose: Calculate the trace of the gas phase rate of strain C
76 ! tensor at i, j, k with cartesian grid modifications C
77 ! C
78 ! Author: Jeff Dietiker Date: 01-Jul-09 C
79 ! C
80 ! C
81 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
82  SUBROUTINE calc_cg_trd_g(IJK,lTRD_G)
83 
84 ! Modules
85 !-----------------------------------------------
86  USE bc
87  USE compar
88  USE cutcell
89  USE fldvar
90  USE functions
91  USE geometry
92  USE indices
93  USE parallel
94  USE param, only: dimension_3
95  USE param1, only: zero, half, undefined
96  USE quadric
97  USE sendrecv
98  IMPLICIT NONE
99 
100 ! Dummy arguments
101 !-----------------------------------------------
102 ! index
103  INTEGER, INTENT(IN) :: IJK
104 ! Strain rate tensor components for mth solids phase
105  DOUBLE PRECISION, INTENT(INOUT) :: ltrD_g(dimension_3)
106 
107 ! Local variables
108 !-----------------------------------------------
109 ! Indices
110  INTEGER :: I, J, K, IMJK, IJMK, IJKM, IM
111 
112  DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
113  DOUBLE PRECISION :: dudx,dvdy,dwdz
114  DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Vi,Wi,Sx,Sy,Sz
115  DOUBLE PRECISION :: UW_g,VW_g,WW_g
116  LOGICAL :: U_NODE_AT_E, U_NODE_AT_W
117  LOGICAL :: V_NODE_AT_N, V_NODE_AT_S
118  LOGICAL :: W_NODE_AT_T, W_NODE_AT_B
119  INTEGER :: BCV
120  INTEGER :: BCT
121 !-----------------------------------------------
122 
123 
124  i = i_of(ijk)
125  j = j_of(ijk)
126  k = k_of(ijk)
127  im = im1(i)
128  imjk = im_of(ijk)
129  ijmk = jm_of(ijk)
130  ijkm = km_of(ijk)
131 
132  IF(flow_at(ijk)) THEN
133  ltrd_g(ijk) = zero
134  RETURN
135  ENDIF
136 
137  bcv = bc_id(ijk)
138  IF(bcv > 0 ) THEN
139  bct = bc_type_enum(bcv)
140  ELSE
141  bct = none
142  ENDIF
143  SELECT CASE (bct)
144  CASE (cg_nsw)
145  noc_trdg = .true.
146  uw_g = zero
147  vw_g = zero
148  ww_g = zero
149  CASE (cg_fsw)
150  noc_trdg = .false.
151  uw_g = zero
152  vw_g = zero
153  ww_g = zero
154  CASE(cg_psw)
155  IF(bc_hw_g(bcv)==undefined) THEN ! same as NSW
156  noc_trdg = .true.
157  uw_g = bc_uw_g(bcv)
158  vw_g = bc_vw_g(bcv)
159  ww_g = bc_ww_g(bcv)
160  ELSEIF(bc_hw_g(bcv)==zero) THEN ! same as FSW
161  noc_trdg = .false.
162  uw_g = zero
163  vw_g = zero
164  ww_g = zero
165  ELSE ! partial slip
166  noc_trdg = .false.
167  ENDIF
168  CASE (cg_mi)
169  ltrd_g(ijk) = zero
170  RETURN
171  CASE (cg_po)
172  ltrd_g(ijk) = zero
173  RETURN
174  CASE (none)
175  ltrd_g(ijk) = zero
176  RETURN
177  END SELECT
178 
179 
180 ! du/dx
181 !=======================================================================
182  u_node_at_e = ((.NOT.blocked_u_cell_at(ijk)) .AND.&
183  (.NOT.wall_u_at(ijk)))
184  u_node_at_w = ((.NOT.blocked_u_cell_at(imjk)).AND.&
185  (.NOT.wall_u_at(imjk)))
186  IF(u_node_at_e.AND.u_node_at_w) THEN
187  ui = half * (u_g(ijk) + u_g(imjk))
188  xi = half * (x_u(ijk) + x_u(imjk))
189  yi = half * (y_u(ijk) + y_u(imjk))
190  zi = half * (z_u(ijk) + z_u(imjk))
191  sx = x_u(ijk) - x_u(imjk)
192  sy = y_u(ijk) - y_u(imjk)
193  sz = z_u(ijk) - z_u(imjk)
194  CALL get_del_h(ijk,'SCALAR',xi,yi,zi,del_h,nx,ny,nz)
195  IF(abs(sx) > zero) THEN
196  dudx = (u_g(ijk) - u_g(imjk))/sx
197  IF(noc_trdg) dudx = dudx - ((ui-uw_g)/(sx*del_h)*&
198  (sy*ny+sz*nz))
199  ELSE
200  dudx = zero
201  ENDIF
202  ELSE
203  dudx = zero
204  ENDIF
205 
206 ! dv/dy
207 !=======================================================================
208  v_node_at_n = ((.NOT.blocked_v_cell_at(ijk)) .AND.&
209  (.NOT.wall_v_at(ijk)))
210  v_node_at_s = ((.NOT.blocked_v_cell_at(ijmk)).AND.&
211  (.NOT.wall_v_at(ijmk)))
212  IF(v_node_at_n.AND.v_node_at_s) THEN
213  vi = half * (v_g(ijk) + v_g(ijmk))
214  xi = half * (x_v(ijk) + x_v(ijmk))
215  yi = half * (y_v(ijk) + y_v(ijmk))
216  zi = half * (z_v(ijk) + z_v(ijmk))
217  sx = x_v(ijk) - x_v(ijmk)
218  sy = y_v(ijk) - y_v(ijmk)
219  sz = z_v(ijk) - z_v(ijmk)
220  CALL get_del_h(ijk,'SCALAR',xi,yi,zi,del_h,nx,ny,nz)
221  IF(abs(sy) > zero) THEN
222  dvdy = (v_g(ijk) - v_g(ijmk))/sy
223  IF(noc_trdg) dvdy = dvdy - ((vi-vw_g)/(sy*del_h)*&
224  (sx*nx+sz*nz))
225  ELSE
226  dvdy = zero
227  ENDIF
228  ELSEIF (v_node_at_n.AND.(.NOT.v_node_at_s).AND.noc_trdg) THEN
229  CALL get_del_h(ijk,'SCALAR',x_v(ijk),y_v(ijk),&
230  z_v(ijk),del_h,nx,ny,nz)
231  dvdy = (v_g(ijk) - vw_g) / del_h * ny
232  ELSEIF ((.NOT.v_node_at_n).AND.v_node_at_s.AND.noc_trdg) THEN
233  CALL get_del_h(ijk,'SCALAR',x_v(ijmk),y_v(ijmk),&
234  z_v(ijmk),del_h,nx,ny,nz)
235  dvdy = (v_g(ijmk) - vw_g) / del_h * ny
236  ELSE
237  dvdy = zero
238  ENDIF
239 
240 ! dw/dz
241 !=======================================================================
242  IF(no_k) THEN
243  dwdz = zero
244  ELSE
245  w_node_at_t = ((.NOT.blocked_w_cell_at(ijk)) .AND.&
246  (.NOT.wall_w_at(ijk)))
247  w_node_at_b = ((.NOT.blocked_w_cell_at(ijkm)).AND.&
248  (.NOT.wall_w_at(ijkm)))
249  IF(w_node_at_t.AND.w_node_at_b) THEN
250  wi = half * (w_g(ijk) + w_g(ijkm))
251  xi = half * (x_w(ijk) + x_w(ijkm))
252  yi = half * (y_w(ijk) + y_w(ijkm))
253  zi = half * (z_w(ijk) + z_w(ijkm))
254  sx = x_w(ijk) - x_w(ijkm)
255  sy = y_w(ijk) - y_w(ijkm)
256  sz = z_w(ijk) - z_w(ijkm)
257  CALL get_del_h(ijk,'SCALAR',xi,yi,zi,del_h,nx,ny,nz)
258  IF(abs(sz) > zero) THEN
259  dwdz = (w_g(ijk) - w_g(ijkm))/sz
260  IF(noc_trdg) dwdz = dwdz - ((wi-ww_g)/(sz*del_h)*&
261  (sx*nx+sy*ny))
262  ELSE
263  dwdz = zero
264  ENDIF
265  ELSEIF (w_node_at_t.AND.(.NOT.w_node_at_b).AND.noc_trdg) THEN
266  CALL get_del_h(ijk,'SCALAR',x_w(ijk),y_w(ijk),&
267  z_w(ijk),del_h,nx,ny,nz)
268  dwdz = (w_g(ijk) - ww_g) / del_h * nz
269  ELSEIF ((.NOT.w_node_at_t).AND.w_node_at_b.AND.noc_trdg) THEN
270  CALL get_del_h(ijk,'SCALAR',x_w(ijkm),y_w(ijkm),&
271  z_w(ijkm),del_h,nx,ny,nz)
272  dwdz = (w_g(ijkm) - ww_g) / del_h * nz
273  ELSE
274  dwdz = zero
275  ENDIF
276  ENDIF ! end if/else no_k
277 
278  ltrd_g(ijk) = dudx + dvdy + dwdz
279 
280  RETURN
281  END SUBROUTINE calc_cg_trd_g
282 
283 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
284 ! C
285 ! Subroutine: GET_FACE_VEL_GAS C
286 ! Purpose: Evaluate the velocity components at each of the faces of C
287 ! a scalar cell C
288 ! C
289 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
290  SUBROUTINE get_face_vel_gas(IJK, uge, ugw, ugn, ugs, ugt, ugb, ugcc, &
291  vge, vgw, vgn, vgs, vgt, vgb, &
292  wge, wgw, wgn, wgs, wgt, wgb, wgcc)
293 
294 ! Modules
295 !-----------------------------------------------------------------------
296  use fldvar, only: u_g, v_g, w_g
297  use functions, only: im_of, jm_of, km_of
298  use functions, only: ip_of, jp_of, kp_of
299  use fun_avg, only: avg_x, avg_x_e
300  use fun_avg, only: avg_y, avg_y_n
301  use fun_avg, only: avg_z, avg_z_t
302  use geometry, only: cylindrical
303  use indices, only: i_of, j_of, k_of
304  use indices, only: im1, jm1, km1
305  use param1, only: zero
306  IMPLICIT NONE
307 
308 ! Dummy arguments
309 !-----------------------------------------------------------------------
310 ! index
311  INTEGER, INTENT(IN) :: IJK
312 ! U_g at the east (i+1/2, j, k) and west face (i-1/2, j, k)
313  DOUBLE PRECISION, INTENT(OUT) :: UgE, UgW
314 ! U_g at the north (i, j+1/2, k) and south face (i, j-1/2, k)
315  DOUBLE PRECISION, INTENT(OUT) :: UgN, UgS
316 ! U_g at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
317  DOUBLE PRECISION, INTENT(OUT) :: UgT, UgB
318 ! U_g at the center of a scalar cell (i, j, k)
319 ! Calculated for Cylindrical coordinates only.
320  DOUBLE PRECISION, INTENT(OUT) :: UgcC
321 
322 ! V_g at the east (i+1/2, j, k) and west face (i-1/2, j, k)
323  DOUBLE PRECISION, INTENT(OUT) :: VgE, VgW
324 ! V_g at the north (i, j+1/2, k) and south face (i, j-1/2, k)
325  DOUBLE PRECISION, INTENT(OUT) :: VgN, VgS
326 ! V_g at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
327  DOUBLE PRECISION, INTENT(OUT) :: VgT, VgB
328 
329 ! W_g at the east (i+1/2, j, k) and west face (i-1/2, j, k)
330  DOUBLE PRECISION, INTENT(OUT) :: WgE, WgW
331 ! W_g at the north (i, j+1/2, k) and south face (i, j-1/2, k)
332  DOUBLE PRECISION, INTENT(OUT) :: WgN, WgS
333 ! W_g at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
334  DOUBLE PRECISION, INTENT(OUT) :: WgT, WgB
335 ! W_g at the center of a scalar cell (i, j, k).
336 ! Calculated for Cylindrical coordinates only.
337  DOUBLE PRECISION, INTENT(OUT) :: WgcC
338 
339 ! Local variables
340 !-----------------------------------------------------------------------
341 ! Cell indices
342  INTEGER :: I, J, K, IM, JM, KM
343  INTEGER :: IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
344  INTEGER :: IMJPK, IMJMK, IMJKP, IMJKM, IPJKM, IPJMK
345  INTEGER :: IJMKP, IJMKM, IJPKM
346 !-----------------------------------------------------------------------
347 
348  i = i_of(ijk)
349  j = j_of(ijk)
350  k = k_of(ijk)
351  im = im1(i)
352  jm = jm1(j)
353  km = km1(k)
354  imjk = im_of(ijk)
355  ipjk = ip_of(ijk)
356  ijmk = jm_of(ijk)
357  ijpk = jp_of(ijk)
358  ijkm = km_of(ijk)
359  ijkp = kp_of(ijk)
360  imjpk = im_of(ijpk)
361  imjmk = im_of(ijmk)
362  imjkp = im_of(ijkp)
363  imjkm = im_of(ijkm)
364  ipjkm = ip_of(ijkm)
365  ipjmk = ip_of(ijmk)
366  ijmkp = jm_of(ijkp)
367  ijmkm = jm_of(ijkm)
368  ijpkm = jp_of(ijkm)
369 
370 ! Find fluid velocity values at faces of the cell
371  ugn = avg_y(avg_x_e(u_g(imjk),u_g(ijk),i), &
372  avg_x_e(u_g(imjpk),u_g(ijpk),i),j) !i, j+1/2, k
373  ugs = avg_y(avg_x_e(u_g(imjmk),u_g(ijmk),i),&
374  avg_x_e(u_g(imjk),u_g(ijk),i),jm) !i, j-1/2, k
375  ugt = avg_z(avg_x_e(u_g(imjk),u_g(ijk),i),&
376  avg_x_e(u_g(imjkp),u_g(ijkp),i),k) !i, j, k+1/2
377  ugb = avg_z(avg_x_e(u_g(imjkm),u_g(ijkm),i),&
378  avg_x_e(u_g(imjk),u_g(ijk),i),km) !i, j, k-1/2
379  uge = u_g(ijk)
380  ugw = u_g(imjk)
381 
382  vge = avg_x(avg_y_n(v_g(ijmk),v_g(ijk)),&
383  avg_y_n(v_g(ipjmk),v_g(ipjk)),i) !i+1/2, j, k
384  vgw = avg_x(avg_y_n(v_g(imjmk),v_g(imjk)),&
385  avg_y_n(v_g(ijmk),v_g(ijk)),im) !i-1/2, j, k
386  vgt = avg_z(avg_y_n(v_g(ijmk),v_g(ijk)),&
387  avg_y_n(v_g(ijmkp),v_g(ijkp)),k) !i, j, k+1/2
388  vgb = avg_z(avg_y_n(v_g(ijmkm),v_g(ijkm)),&
389  avg_y_n(v_g(ijmk),v_g(ijk)),km) !i, j, k-1/2
390  vgn = v_g(ijk)
391  vgs = v_g(ijmk)
392 
393  wgn = avg_y(avg_z_t(w_g(ijkm),w_g(ijk)),&
394  avg_z_t(w_g(ijpkm),w_g(ijpk)),j) !i, j+1/2, k
395  wgs = avg_y(avg_z_t(w_g(ijmkm),w_g(ijmk)),&
396  avg_z_t(w_g(ijkm),w_g(ijk)),jm) !i, j-1/2, k
397  wge = avg_x(avg_z_t(w_g(ijkm),w_g(ijk)),&
398  avg_z_t(w_g(ipjkm),w_g(ipjk)),i) !i+1/2, j, k
399  wgw = avg_x(avg_z_t(w_g(imjkm),w_g(imjk)),&
400  avg_z_t(w_g(ijkm),w_g(ijk)),im) !i-1/2, j, k
401  wgt = w_g(ijk)
402  wgb = w_g(ijkm)
403 
404  IF (cylindrical) THEN
405  ugcc = avg_x_e(u_g(imjk),u_g(ijk),i) !i, j, k
406  wgcc = avg_z_t(w_g(ijkm),w_g(ijk)) !i, j, k
407  ELSE
408  ugcc = zero
409  wgcc = zero
410  ENDIF
411 
412  RETURN
413  END SUBROUTINE get_face_vel_gas
414 
415 
416 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
417 ! C
418 ! Subroutine: CALC_DERIV_VEL_GAS C
419 ! Purpose: Calculate the gradient of the gas phase velocity and C
420 ! the rate of strain tensor at i, j, k C
421 ! C
422 ! | du/dx du/dy du/dz | C
423 ! DELV = | dv/dx dv/dy dv/dz | = dVi/dxj C
424 ! | dw/dx dw/dy dw/dz | C
425 ! C
426 ! 1/2(DELV+DELV^T) = 1/2(dVi/dxj+dVj/dxi) C
427 ! C
428 ! C
429 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
430  SUBROUTINE calc_deriv_vel_gas(IJK, lVelGradG, lRateStrainG)
431 
432 ! Modules
433 !-----------------------------------------------------------------------
434  use cutcell, only: cut_cell_at
435  use geometry, only: odx, ody, odz
436  use geometry, only: ox
437  use indices, only: i_of, j_of, k_of
438  use param1, only: half
439 ! Dummy arguments
440 !-----------------------------------------------------------------------
441 ! index
442  INTEGER, INTENT(IN) :: IJK
443 ! gradient of velocity
444  DOUBLE PRECISION, INTENT(OUT) :: lVelGradG(3,3) ! delV
445 ! rate of strain tensor
446  DOUBLE PRECISION, INTENT(OUT) :: lRateStrainG(3,3) ! D_g
447 
448 ! Local variables
449 !-----------------------------------------------------------------------
450 ! Cell indices
451  INTEGER :: I, J, K
452 ! face values of velocity
453  DOUBLE PRECISION :: uge, ugw, ugn, ugs, ugt, ugb, ugcc
454  DOUBLE PRECISION :: vge, vgw, vgn, vgs, vgt, vgb
455  DOUBLE PRECISION :: wge, wgw, wgn, wgs, wgt, wgb, wgcc
456 !-----------------------------------------------------------------------
457 
458  IF(.NOT.cut_cell_at(ijk)) THEN
459  i = i_of(ijk)
460  j = j_of(ijk)
461  k = k_of(ijk)
462 
463 ! Set the velocity components at the scalar faces
464  CALL get_face_vel_gas(ijk, uge, ugw, ugn, ugs, ugt, ugb, ugcc, &
465  vge, vgw, vgn, vgs, vgt, vgb, &
466  wge, wgw, wgn, wgs, wgt, wgb, wgcc)
467 
468 ! Gradient of gas velocity at cell center (i, j, k)
469  lvelgradg(1,1) = (uge-ugw)*odx(i)
470  lvelgradg(1,2) = (ugn-ugs)*ody(j)
471  lvelgradg(1,3) = (ugt-ugb)*(ox(i)*odz(k))-wgcc*ox(i)
472 
473  lvelgradg(2,1) = (vge-vgw)*odx(i)
474  lvelgradg(2,2) = (vgn-vgs)*ody(j)
475  lvelgradg(2,3) = (vgt-vgb)*(ox(i)*odz(k))
476 
477  lvelgradg(3,1) = (wge-wgw)*odx(i)
478  lvelgradg(3,2) = (wgn-wgs)*ody(j)
479  lvelgradg(3,3) = (wgt-wgb)*(ox(i)*odz(k)) + ugcc*ox(i)
480 
481 
482 ! Strain rate tensor, D_G, at cell center (i, j, k)
483 ! or calculated as 0.5*(DelG+DelG^T)
484  lratestraing(1,1) = (uge-ugw)*odx(i)
485  lratestraing(1,2) = half*((ugn-ugs)*ody(j)+&
486  (vge-vgw)*odx(i))
487  lratestraing(1,3) = half*((wge-wgw)*odx(i)+&
488  (ugt-ugb)*(ox(i)*odz(k))-&
489  wgcc*ox(i))
490 
491  lratestraing(2,1) = lratestraing(1,2)
492  lratestraing(2,2) = (vgn-vgs)*ody(j)
493  lratestraing(2,3) = half*((vgt-vgb)*(ox(i)*odz(k))+&
494  (wgn-wgs)*ody(j))
495 
496  lratestraing(3,1) = lratestraing(1,3)
497  lratestraing(3,2) = lratestraing(2,3)
498  lratestraing(3,3) = (wgt-wgb)*(ox(i)*odz(k)) +&
499  ugcc*ox(i)
500 
501  ELSE ! cut cell
502  CALL calc_cg_deriv_vel_gas(ijk,lvelgradg,lratestraing)
503  ENDIF
504 
505  RETURN
506  END SUBROUTINE calc_deriv_vel_gas
507 
508 
509 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
510 ! C
511 ! Module name: CALC_CG_DERIV_VEL_GAS C
512 ! Purpose: Calculate velocity derivatives in scalar cut-cell C
513 ! gas phase C
514 ! C
515 ! Author: Jeff Dietiker Date: 25-JAN-96 C
516 ! C
517 ! C
518 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
519  SUBROUTINE calc_cg_deriv_vel_gas(IJK, lVelGradG, lRateStrainG)
520 
521 ! Modules
522 !-----------------------------------------------
523  USE param
524  USE param1
525  USE parallel
526  USE geometry
527  USE fldvar
528  USE indices
529  USE compar
530  USE sendrecv
531  USE bc
532  USE cutcell
533  USE quadric
534  USE functions
535  IMPLICIT NONE
536 
537 ! Dummy arguments
538 !-----------------------------------------------
539 ! index
540  INTEGER, INTENT(IN) :: IJK
541 ! velocity gradient gas phase
542 ! | du/dx du/dy du/dz |
543 ! DELV = | dv/dx dv/dy dv/dz | = dUi/dxj
544 ! | dw/dx dw/dy dw/dz |
545  DOUBLE PRECISION, INTENT(OUT) :: lVelGradG(3,3)
546 ! rate of strain tensor gas phase
547  DOUBLE PRECISION, INTENT(OUT) :: lRateStrainG(3,3)
548 
549 ! Local variables
550 !-----------------------------------------------
551 ! Indices
552  INTEGER :: I, J, K, IMJK, IJMK, IJKM
553 ! loop counters
554  INTEGER :: P, Q
555 !
556  DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
557  DOUBLE PRECISION :: dudx,dudy,dudz
558  DOUBLE PRECISION :: dvdx,dvdy,dvdz
559  DOUBLE PRECISION :: dwdx,dwdy,dwdz
560  DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Vi,Wi,Sx,Sy,Sz
561  DOUBLE PRECISION :: UW_g,VW_g,WW_g
562 
563  LOGICAL :: U_NODE_AT_E, U_NODE_AT_W
564  LOGICAL :: V_NODE_AT_N, V_NODE_AT_S
565  LOGICAL :: W_NODE_AT_T, W_NODE_AT_B
566  INTEGER :: BCV
567  INTEGER :: BCT
568 !-----------------------------------------------
569 
570 ! initialize
571  lvelgradg = zero
572  lratestraing = zero
573 
574  i = i_of(ijk)
575  j = j_of(ijk)
576  k = k_of(ijk)
577  imjk = im_of(ijk)
578  ijmk = jm_of(ijk)
579  ijkm = km_of(ijk)
580 
581  IF(flow_at(ijk)) THEN
582  lvelgradg = zero
583  RETURN
584  ENDIF
585 
586  bcv = bc_id(ijk)
587  IF(bcv > 0 ) THEN
588  bct = bc_type_enum(bcv)
589  ELSE
590  bct = none
591  ENDIF
592  SELECT CASE (bct)
593  CASE (cg_nsw)
594  noc_trdg = .true.
595  uw_g = zero
596  vw_g = zero
597  ww_g = zero
598  CASE (cg_fsw)
599  noc_trdg = .false.
600  uw_g = zero
601  vw_g = zero
602  ww_g = zero
603  CASE(cg_psw)
604  IF(bc_hw_g(bcv)==undefined) THEN ! same as NSW
605  noc_trdg = .true.
606  uw_g = bc_uw_g(bcv)
607  vw_g = bc_vw_g(bcv)
608  ww_g = bc_ww_g(bcv)
609  ELSEIF(bc_hw_g(bcv)==zero) THEN ! same as FSW
610  noc_trdg = .false.
611  uw_g = zero
612  vw_g = zero
613  ww_g = zero
614  ELSE ! partial slip
615  noc_trdg = .false.
616  ENDIF
617  CASE (cg_mi)
618  lvelgradg = zero
619  RETURN
620  CASE (cg_po)
621  lvelgradg = zero
622  RETURN
623  CASE (none)
624  lvelgradg = zero
625  RETURN
626  END SELECT
627 
628 
629 ! du/dx, du/dy, du/dz
630 !=======================================================================
631  u_node_at_e = ((.NOT.blocked_u_cell_at(ijk)) .AND.&
632  (.NOT.wall_u_at(ijk)))
633  u_node_at_w = ((.NOT.blocked_u_cell_at(imjk)).AND.&
634  (.NOT.wall_u_at(imjk)))
635 
636  IF(u_node_at_e.AND.u_node_at_w) THEN
637  ui = half * (u_g(ijk) + u_g(imjk))
638  xi = half * (x_u(ijk) + x_u(imjk))
639  yi = half * (y_u(ijk) + y_u(imjk))
640  zi = half * (z_u(ijk) + z_u(imjk))
641  sx = x_u(ijk) - x_u(imjk)
642  sy = y_u(ijk) - y_u(imjk)
643  sz = z_u(ijk) - z_u(imjk)
644  CALL get_del_h(ijk,'SCALAR',xi,yi,zi,del_h,nx,ny,nz)
645  IF(abs(sx) > zero) THEN
646  dudx = (u_g(ijk) - u_g(imjk))/sx
647  dudy = zero
648  dudz = zero
649  IF(noc_trdg) THEN
650  dudx = dudx - ((ui-uw_g)/(sx*del_h)*(sy*ny+sz*nz))
651  dudy = (ui-uw_g) / del_h * ny
652  dudz = (ui-uw_g) / del_h * nz
653  ENDIF
654  ELSE
655  dudx = zero
656  dudy = zero
657  dudz = zero
658  ENDIF
659  ELSEIF (u_node_at_e.AND.(.NOT.u_node_at_w).AND.noc_trdg) THEN
660  CALL get_del_h(ijk,'SCALAR',x_u(ijk),y_u(ijk),z_u(ijk),&
661  del_h,nx,ny,nz)
662  dudx = (u_g(ijk) - uw_g) / del_h * nx
663  dudy = (u_g(ijk) - uw_g) / del_h * ny
664  dudz = (u_g(ijk) - uw_g) / del_h * nz
665  ELSEIF ((.NOT.u_node_at_e).AND.u_node_at_w.AND.noc_trdg) THEN
666  CALL get_del_h(ijk,'SCALAR',x_u(imjk),y_u(imjk),z_u(imjk),&
667  del_h,nx,ny,nz)
668  dudx = (u_g(imjk) - uw_g) / del_h * nx
669  dudy = (u_g(imjk) - uw_g) / del_h * ny
670  dudz = (u_g(imjk) - uw_g) / del_h * nz
671  ELSE
672  dudx = zero
673  dudy = zero
674  dudz = zero
675  ENDIF
676  lvelgradg(1,1) = dudx
677  lvelgradg(1,2) = dudy
678  lvelgradg(1,3) = dudz
679 
680 ! dv/dx, dv/dy, dv/dz
681 !=======================================================================
682  v_node_at_n = ((.NOT.blocked_v_cell_at(ijk)) .AND.&
683  (.NOT.wall_v_at(ijk)))
684  v_node_at_s = ((.NOT.blocked_v_cell_at(ijmk)).AND.&
685  (.NOT.wall_v_at(ijmk)))
686 
687  IF(v_node_at_n.AND.v_node_at_s) THEN
688  vi = half * (v_g(ijk) + v_g(ijmk))
689  xi = half * (x_v(ijk) + x_v(ijmk))
690  yi = half * (y_v(ijk) + y_v(ijmk))
691  zi = half * (z_v(ijk) + z_v(ijmk))
692  sx = x_v(ijk) - x_v(ijmk)
693  sy = y_v(ijk) - y_v(ijmk)
694  sz = z_v(ijk) - z_v(ijmk)
695  CALL get_del_h(ijk,'SCALAR',xi,yi,zi,del_h,nx,ny,nz)
696  IF(abs(sy) > zero) THEN
697  dvdx = zero
698  dvdy = (v_g(ijk) - v_g(ijmk))/sy
699  dvdz = zero
700  IF(noc_trdg) THEN
701  dvdx = (vi-vw_g) / del_h * nx
702  dvdy = dvdy - ((vi-vw_g)/(sy*del_h)*(sx*nx+sz*nz))
703  dvdz = (vi-vw_g) / del_h * nz
704  ENDIF
705  ELSE
706  dvdx = zero
707  dvdy = zero
708  dvdz = zero
709  ENDIF
710  ELSEIF (v_node_at_n.AND.(.NOT.v_node_at_s).AND.noc_trdg) THEN
711  CALL get_del_h(ijk,'SCALAR',x_v(ijk),y_v(ijk),z_v(ijk),&
712  del_h,nx,ny,nz)
713  dvdx = (v_g(ijk) - vw_g) / del_h * nx
714  dvdy = (v_g(ijk) - vw_g) / del_h * ny
715  dvdz = (v_g(ijk) - vw_g) / del_h * nz
716  ELSEIF ((.NOT.v_node_at_n).AND.v_node_at_s.AND.noc_trdg) THEN
717  CALL get_del_h(ijk,'SCALAR',x_v(ijmk),y_v(ijmk),z_v(ijmk),&
718  del_h,nx,ny,nz)
719  dvdx = (v_g(ijmk) - vw_g) / del_h * nx
720  dvdy = (v_g(ijmk) - vw_g) / del_h * ny
721  dvdz = (v_g(ijmk) - vw_g) / del_h * nz
722  ELSE
723  dvdx = zero
724  dvdy = zero
725  dvdz = zero
726  ENDIF
727  lvelgradg(2,1) = dvdx
728  lvelgradg(2,2) = dvdy
729  lvelgradg(2,3) = dvdz
730 
731 ! dw/dx, dw/dy, dw/dz
732 !=======================================================================
733  IF(no_k) THEN
734  dwdx = zero
735  dwdy = zero
736  dwdz = zero
737  ELSE
738  w_node_at_t = ((.NOT.blocked_w_cell_at(ijk)) .AND.&
739  (.NOT.wall_w_at(ijk)))
740  w_node_at_b = ((.NOT.blocked_w_cell_at(ijkm)).AND.&
741  (.NOT.wall_w_at(ijkm)))
742  IF(w_node_at_t.AND.w_node_at_b) THEN
743  wi = half * (w_g(ijk) + w_g(ijkm))
744  xi = half * (x_w(ijk) + x_w(ijkm))
745  yi = half * (y_w(ijk) + y_w(ijkm))
746  zi = half * (z_w(ijk) + z_w(ijkm))
747  sx = x_w(ijk) - x_w(ijkm)
748  sy = y_w(ijk) - y_w(ijkm)
749  sz = z_w(ijk) - z_w(ijkm)
750  CALL get_del_h(ijk,'SCALAR',xi,yi,zi,del_h,nx,ny,nz)
751  IF(abs(sz) > zero) THEN
752  dwdx = zero
753  dwdy = zero
754  dwdz = (w_g(ijk) - w_g(ijkm))/sz
755  IF(noc_trdg) THEN
756  dwdx = (wi-ww_g) / del_h * nx
757  dwdy = (wi-ww_g) / del_h * ny
758  dwdz = dwdz - ((wi-ww_g)/(sz*del_h)*(sx*nx+sy*ny))
759  ENDIF
760  ELSE
761  dwdx = zero
762  dwdy = zero
763  dwdz = zero
764  ENDIF
765  ELSEIF (w_node_at_t.AND.(.NOT.w_node_at_b).AND.noc_trdg) THEN
766  CALL get_del_h(ijk,'SCALAR',x_w(ijk),y_w(ijk),z_w(ijk),&
767  del_h,nx,ny,nz)
768  dwdx = (w_g(ijk) - ww_g) / del_h * nx
769  dwdy = (w_g(ijk) - ww_g) / del_h * ny
770  dwdz = (w_g(ijk) - ww_g) / del_h * nz
771  ELSEIF ((.NOT.w_node_at_t).AND.w_node_at_b.AND.noc_trdg) THEN
772  CALL get_del_h(ijk,'SCALAR',x_w(ijkm),y_w(ijkm),z_w(ijkm),&
773  del_h,nx,ny,nz)
774  dwdx = (w_g(ijkm) - ww_g) / del_h * nx
775  dwdy = (w_g(ijkm) - ww_g) / del_h * ny
776  dwdz = (w_g(ijkm) - ww_g) / del_h * nz
777  ELSE
778  dwdx = zero
779  dwdy = zero
780  dwdz = zero
781  ENDIF
782  ENDIF
783  lvelgradg(3,1) = dwdx
784  lvelgradg(3,2) = dwdy
785  lvelgradg(3,3) = dwdz
786 
787  DO p = 1,3
788  DO q = 1,3
789  lratestraing(p,q) = half * (lvelgradg(p,q)+lvelgradg(q,p))
790  ENDDO
791  ENDDO
792 
793 
794  RETURN
795  END SUBROUTINE calc_cg_deriv_vel_gas
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
geometry::ody
double precision, dimension(:), allocatable ody
Definition: geometry_mod.f:116
geometry::odx
double precision, dimension(:), allocatable odx
Definition: geometry_mod.f:114
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
get_face_vel_gas
subroutine get_face_vel_gas(IJK, uge, ugw, ugn, ugs, ugt, ugb, ugc
Definition: calc_trd_g.f:291
geometry::x_e
double precision, dimension(:), allocatable x_e
Definition: geometry_mod.f:134
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::wall_v_at
logical, dimension(:), allocatable wall_v_at
Definition: cutcell_mod.f:127
calc_cg_deriv_vel_gas
subroutine calc_cg_deriv_vel_gas(IJK, lVelGradG, lRateStrainG)
Definition: calc_trd_g.f:520
bc::bc_type_enum
integer, dimension(dimension_bc) bc_type_enum
Definition: bc_mod.f:146
calc_trd_g
subroutine calc_trd_g(lTRD_G)
Definition: calc_trd_g.f:11
param1::undefined
double precision, parameter undefined
Definition: param1_mod.f:18
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
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
calc_deriv_vel_gas
subroutine calc_deriv_vel_gas(IJK, lVelGradG, lRateStrainG)
Definition: calc_trd_g.f:431
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
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
geometry::ox
double precision, dimension(:), allocatable ox
Definition: geometry_mod.f:140
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
fldvar::v_g
double precision, dimension(:), allocatable v_g
Definition: fldvar_mod.f:99
cutcell::noc_trdg
logical noc_trdg
Definition: cutcell_mod.f:389
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
param
Definition: param_mod.f:2
geometry::no_k
logical no_k
Definition: geometry_mod.f:28
cutcell::cut_cell_at
logical, dimension(:), allocatable cut_cell_at
Definition: cutcell_mod.f:355
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
cutcell::bc_id
integer, dimension(:), allocatable bc_id
Definition: cutcell_mod.f:433
fldvar::u_g
double precision, dimension(:), allocatable u_g
Definition: fldvar_mod.f:87
cutcell::blocked_v_cell_at
logical, dimension(:), allocatable blocked_v_cell_at
Definition: cutcell_mod.f:365
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
bc::bc_ww_g
double precision, dimension(dimension_bc) bc_ww_g
Definition: bc_mod.f:319
calc_cg_trd_g
subroutine calc_cg_trd_g(IJK, lTRD_G)
Definition: calc_trd_g.f:83
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