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