MFIX  2016-1
calc_grbdry.f
Go to the documentation of this file.
1 MODULE calc_gr_boundary
2  USE exit, only: mfix_exit
3  CONTAINS
4 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
5 ! C
6 ! Subroutine: CALC_GRBDRY C
7 ! Purpose: Main controller subroutine for calculating coefficients C
8 ! for momentum boundary conditions using kinetic & frictional C
9 ! theory C
10 ! C
11 ! Author: K. Agrawal & A. Srivastava, Princeton Univ. Date: 19-JAN-98 C
12 ! Reviewer: Date: C
13 ! C
14 ! C
15 ! Literature/Document References: C
16 ! C
17 ! C
18 ! Notes: C
19 ! This routine will not be entered unless granular_energy so we C
20 ! do not need to check for its condition! C
21 ! C
22 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
23 
24  SUBROUTINE calc_grbdry(IJK1, IJK2, FCELL, COM, M, L, Gw, Hw, Cw)
25 
26 !-----------------------------------------------
27 ! Modules
28 !-----------------------------------------------
29  USE bc
30  USE compar
31  USE constant
32  USE fldvar
33  USE fun_avg
34  USE functions
35  USE geometry
36  USE indices
37  USE mpi_utility
38  USE param
39  USE param1
40  USE physprop
41  USE rdf
42  USE run
43  USE toleranc
44  USE turb
45  USE visc_s
46  IMPLICIT NONE
47 
48 !-----------------------------------------------
49 ! Dummy Arguments
50 !-----------------------------------------------
51 ! IJK indices for wall cell (ijk1) and fluid cell (ijk2)
52  INTEGER, INTENT(IN) :: IJK1, IJK2
53 ! The location (e,w,n...) of fluid cell
54  CHARACTER, INTENT(IN) :: FCELL
55 ! Velocity component (U, V, W)
56  CHARACTER :: COM
57 ! Solids phase index
58  INTEGER, INTENT(IN) :: M
59 ! Index corresponding to boundary condition
60  INTEGER, INTENT(IN) :: L
61 ! Wall momentum coefficients:
62 ! 1st term on LHS
63  DOUBLE PRECISION, INTENT(INOUT) :: Gw
64 ! 2nd term on LHS
65  DOUBLE PRECISION, INTENT(INOUT) :: Hw
66 ! all terms appearing on RHS
67  DOUBLE PRECISION, INTENT(INOUT) :: Cw
68 !-----------------------------------------------
69 ! Local Variables
70 !-----------------------------------------------
71 ! Other indices
72  INTEGER :: IJK2E, IPJK2, IPJKM2, IPJKP2, IPJMK2, IPJPK2
73  INTEGER :: IJK2N, IJPK2, IJPKM2, IJPKP2, IMJPK2
74  INTEGER :: IJK2T, IJKP2, IJMKP2, IMJKP2
75 ! Solids phase index
76  INTEGER :: MM
77 !
78  DOUBLE PRECISION :: smallTheta
79 ! Average scalars
80  DOUBLE PRECISION :: EPg_avg, Mu_g_avg, RO_g_avg
81 ! Average void fraction at packing
82  DOUBLE PRECISION :: ep_star_avg
83 ! Average scalars modified to include all solid phases
84  DOUBLE PRECISION :: EPs_avg(dimension_m), DP_avg(dimension_m),&
85  TH_avg(dimension_m), AVGX1, AVGX2
86  DOUBLE PRECISION ROS_AVG(dimension_m)
87 ! Average Simonin and Ahmadi variables (sof)
88  DOUBLE PRECISION :: K_12_avg, Tau_12_avg, Tau_1_avg
89 ! Average velocities
90 ! values of U_sm, V_sm, W_sm at appropriate place on boundary wall
91  DOUBLE PRECISION :: USCM, VSCM,WSCM
92  DOUBLE PRECISION :: USCM1,USCM2,VSCM1,VSCM2,WSCM1,WSCM2
93 ! values of U_g, V_g, W_g at appropriate place on boundary wall
94  DOUBLE PRECISION :: UGC, VGC, WGC
95  DOUBLE PRECISION :: WGC1, WGC2, VGC1, VGC2, UGC1, UGC2
96 ! velocity variables used to standarize dummy argument for different dirs
97  DOUBLE PRECISION :: WVELS
98  DOUBLE PRECISION :: VELS
99 ! del.u
100  DOUBLE PRECISION :: DEL_DOT_U
101 ! S:S
102  DOUBLE PRECISION :: S_DDOT_S
103 ! S_dd (d can be x, y or z)
104  DOUBLE PRECISION :: S_dd
105 ! Magnitude of gas-solids relative velocity
106  DOUBLE PRECISION :: VREL
107 ! slip velocity between wall and particles for Jenkins bc (sof)
108  DOUBLE PRECISION :: VSLIP
109 ! radial distribution function at contact
110  DOUBLE PRECISION :: g0(dimension_m)
111 ! Sum of eps*G_0
112  DOUBLE PRECISION :: g0EPs_avg
113 ! Error message
114  CHARACTER(LEN=80) :: LINE
115 
116 !-----------------------------------------------
117 ! Note: EP_s, MU_g, and RO_g are undefined at IJK1 (wall cell).
118 ! Hence IJK2 (fluid cell) is used in averages.
119 
120  smalltheta = (to_si)**4 * zero_ep_s
121 
122 !-----------------------------------------------
123 ! Calculations for U momentum equation
124  IF (com .EQ. 'U')THEN
125 
126  ijk2e = east_of(ijk2)
127  ipjk2 = ip_of(ijk2)
128 
129  epg_avg = avg_x(ep_g(ijk2), ep_g(ijk2e), i_of(ijk2))
130  ep_star_avg = avg_x(ep_star_array(ijk2), ep_star_array(ijk2e), i_of(ijk2))
131  mu_g_avg = avg_x(mu_g(ijk2), mu_g(ijk2e), i_of(ijk2))
132  ro_g_avg = avg_x(ro_g(ijk2), ro_g(ijk2e), i_of(ijk2))
133 
134  k_12_avg = zero
135  tau_12_avg = zero
136  tau_1_avg = zero
137  IF(kt_type_enum == simonin_1996 .OR. &
138  kt_type_enum == ahmadi_1995) THEN
139  tau_1_avg = avg_x(tau_1(ijk2), tau_1(ijk2e), i_of(ijk2))
140 ! Simonin only:
141  k_12_avg = avg_x(k_12(ijk2), k_12(ijk2e), i_of(ijk2))
142  tau_12_avg = avg_x(tau_12(ijk2), tau_12(ijk2e), i_of(ijk2))
143  ENDIF
144 
145  g0eps_avg = zero
146  DO mm = 1, smax
147  g0(mm) = g_0avg(ijk2, ijk2e, 'X', i_of(ijk2), m, mm)
148  eps_avg(mm) = avg_x(ep_s(ijk2, mm), ep_s(ijk2e, mm), i_of(ijk2))
149  dp_avg(mm) = avg_x(d_p(ijk2,mm), d_p(ijk2e,mm), i_of(ijk2))
150  g0eps_avg = g0eps_avg + g_0avg(ijk2, ijk2e, 'X', i_of(ijk2), m, mm) &
151  * avg_x(ep_s(ijk2, mm), ep_s(ijk2e, mm), i_of(ijk2))
152  ros_avg(mm) = avg_x(ro_s(ijk2, mm), ro_s(ijk2e, mm), i_of(ijk2))
153 
154  th_avg(mm) = avg_x(theta_m(ijk2,mm), theta_m(ijk2e,mm), i_of(ijk2))
155  IF(th_avg(mm) < zero) th_avg(mm) = smalltheta
156  ENDDO
157 
158  wvels = bc_uw_s(l,m)
159 
160 
161  IF(fcell .EQ. 'N')THEN
162  ipjmk2 = jm_of(ipjk2)
163 
164 ! code modified for some corner cells
165  DO mm = 1, smax
166  avgx1 = avg_x(theta_m(ijk1,mm), theta_m(ipjmk2,mm), i_of(ijk1))
167  avgx2 = avg_x(theta_m(ijk2,mm), theta_m(ipjk2,mm), i_of(ijk2))
168  IF(avgx1 < zero .AND. avgx2 > zero) avgx1 = avgx2
169  IF(avgx2 < zero .AND. avgx1 > zero) avgx2 = avgx1
170  IF(avgx1 < zero .AND. avgx2 < zero) THEN
171  th_avg(mm) = smalltheta
172  ELSE
173  th_avg(mm) = avg_y(avgx1, avgx2, j_of(ijk1))
174  ENDIF
175  ENDDO
176 
177 ! Calculate velocity components at i+1/2, j+1/2, k (relative to IJK1)
178  ugc = avg_y(u_g(ijk1), u_g(ijk2),j_of(ijk1))
179  vgc = avg_x(v_g(ijk1), v_g(ipjmk2),i_of(ijk1))
180  wgc1 = avg_x(avg_z_t(w_g(km_of(ijk2)), w_g(ijk2)),&
181  avg_z_t(w_g(km_of(ipjk2)), w_g(ipjk2)),&
182  i_of(ijk2))
183  wgc2 = avg_x(avg_z_t(w_g(km_of(ijk1)), w_g(ijk1)),&
184  avg_z_t(w_g(km_of(ipjmk2)), w_g(ipjmk2)),&
185  i_of(ijk1))
186  wgc = avg_y(wgc2, wgc1, j_of(ijk1))
187  uscm = avg_y(u_s(ijk1,m), u_s(ijk2,m),j_of(ijk1))
188  vscm = avg_x(v_s(ijk1,m), v_s(ipjmk2,m),i_of(ijk1))
189  wscm1= avg_x(avg_z_t(w_s(km_of(ijk2),m),w_s(ijk2,m)),&
190  avg_z_t(w_s(km_of(ipjk2),m),w_s(ipjk2,m)),&
191  i_of(ijk2))
192  wscm2= avg_x(avg_z_t(w_s(km_of(ijk1),m),w_s(ijk1,m)),&
193  avg_z_t(w_s(km_of(ipjmk2),m),w_s(ipjmk2,m)),&
194  i_of(ijk1))
195  wscm = avg_y(wscm2, wscm1, j_of(ijk1))
196  vels = uscm
197 
198  ELSEIF(fcell .EQ. 'S')THEN
199  ipjpk2= jp_of(ipjk2)
200 
201 ! code modified for some corner cells
202  DO mm = 1, smax
203  avgx1 = avg_x(theta_m(ijk2,mm),theta_m(ipjk2,mm),i_of(ijk2))
204  avgx2 = avg_x(theta_m(ijk1,mm),theta_m(ipjpk2,mm),i_of(ijk1))
205  IF(avgx1 < zero .AND. avgx2 > zero) avgx1 = avgx2
206  IF(avgx2 < zero .AND. avgx1 > zero) avgx2 = avgx1
207  IF(avgx1 < zero .AND. avgx2 < zero) THEN
208  th_avg(mm) = smalltheta
209  ELSE
210  th_avg(mm) = avg_y(avgx1, avgx2, j_of(ijk2))
211  ENDIF
212  ENDDO
213 
214 
215 ! Calculate velocity components at i+1/2, j+1/2, k relative to IJK2
216  ugc = avg_y(u_g(ijk2),u_g(ijk1),j_of(ijk2))
217  vgc = avg_x(v_g(ijk2),v_g(ipjk2),i_of(ijk2))
218  wgc1 = avg_x(avg_z_t(w_g(km_of(ijk2)), w_g(ijk2)),&
219  avg_z_t(w_g(km_of(ipjk2)), w_g(ipjk2)),&
220  i_of(ijk2))
221  wgc2 = avg_x(avg_z_t(w_g(km_of(ijk1)), w_g(ijk1)),&
222  avg_z_t(w_g(km_of(ipjpk2)), w_g(ipjpk2)),&
223  i_of(ijk1))
224  wgc = avg_y(wgc1, wgc2, j_of(ijk2))
225  uscm = avg_y(u_s(ijk2, m),u_s(ijk1, m),j_of(ijk2))
226  vscm = avg_x(v_s(ijk2, m),v_s(ipjk2, m),i_of(ijk2))
227  wscm1= avg_x(avg_z_t(w_s(km_of(ijk2),m),w_s(ijk2,m)),&
228  avg_z_t(w_s(km_of(ipjk2),m),w_s(ipjk2,m)),&
229  i_of(ijk2))
230  wscm2= avg_x(avg_z_t(w_s(km_of(ijk1),m),w_s(ijk1,m)),&
231  avg_z_t(w_s(km_of(ipjpk2),m),w_s(ipjpk2,m)),&
232  i_of(ijk1))
233  wscm = avg_y(wscm1, wscm2, j_of(ijk2))
234  vels = uscm
235 
236  ELSEIF(fcell .EQ. 'T')THEN
237  ipjkm2= km_of(ipjk2)
238 
239 ! code modified for some corner cells
240  DO mm = 1, smax
241  avgx1 = avg_x(theta_m(ijk1,mm),theta_m(ipjkm2,mm),i_of(ijk1))
242  avgx2 = avg_x(theta_m(ijk2,mm),theta_m(ipjk2,mm),i_of(ijk2))
243  IF(avgx1 < zero .AND. avgx2 > zero) avgx1 = avgx2
244  IF(avgx2 < zero .AND. avgx1 > zero) avgx2 = avgx1
245  IF(avgx1 < zero .AND. avgx2 < zero) THEN
246  th_avg(mm) = smalltheta
247  ELSE
248  th_avg(mm) = avg_z(avgx1, avgx2, k_of(ijk1))
249  ENDIF
250  ENDDO
251 
252 ! Calculate velocity components at i+1/2,j,k-1/2 relative to IJK2
253  ugc = avg_z(u_g(ijk1), u_g(ijk2), k_of(ijk1))
254  vgc1 = avg_x(avg_y_n(v_g(jm_of(ijk2)),v_g(ijk2)),&
255  avg_y_n(v_g(jm_of(ipjk2)),v_g(ipjk2)),&
256  i_of(ijk2))
257  vgc2 = avg_x(avg_y_n(v_g(jm_of(ijk1)),v_g(ijk1)),&
258  avg_y_n(v_g(jm_of(ipjkm2)),v_g(ipjkm2)),&
259  i_of(ijk1))
260  vgc = avg_z(vgc2,vgc1,k_of(ijk1))
261  wgc = avg_x(w_g(ijk1), w_g(ipjkm2),i_of(ijk1))
262  uscm = avg_z(u_s(ijk1,m), u_s(ijk2,m), k_of(ijk1))
263  vscm1= avg_x(avg_y_n(v_s(jm_of(ijk2),m),v_s(ijk2,m)),&
264  avg_y_n(v_s(jm_of(ipjk2),m),v_s(ipjk2,m)),&
265  i_of(ijk2))
266  vscm2= avg_x(avg_y_n(v_s(jm_of(ijk1),m),v_s(ijk1,m)),&
267  avg_y_n(v_s(jm_of(ipjkm2),m),v_s(ipjkm2,m)),&
268  i_of(ijk1))
269  vscm = avg_z(vscm2,vscm1,k_of(ijk1))
270  wscm = avg_x(w_s(ijk1,m), w_s(ipjkm2,m), i_of(ijk1))
271  vels = uscm
272 
273  ELSEIF(fcell .EQ. 'B')THEN
274  ipjkp2= kp_of(ipjk2)
275 
276 ! code modified for some corner cells
277  DO mm = 1, smax
278  avgx1 = avg_x(theta_m(ijk1,mm), theta_m(ipjkp2,mm),i_of(ijk1))
279  avgx2 = avg_x(theta_m(ijk2,mm), theta_m(ipjk2,mm),i_of(ijk2))
280  IF(avgx1 < zero .AND. avgx2 > zero) avgx1 = avgx2
281  IF(avgx2 < zero .AND. avgx1 > zero) avgx2 = avgx1
282  IF(avgx1 < zero .AND. avgx2 < zero) THEN
283  th_avg(mm) = smalltheta
284  ELSE
285  th_avg(mm) = avg_z(avgx1, avgx2, k_of(ijk2))
286  ENDIF
287  ENDDO
288 
289 ! Calculate velocity components at i+1/2,j,k-1/2 relative to IJK1
290  ugc = avg_z(u_g(ijk2), u_g(ijk1), k_of(ijk2))
291  vgc1 = avg_x(avg_y_n(v_g(jm_of(ijk2)),v_g(ijk2)),&
292  avg_y_n(v_g(jm_of(ipjk2)),v_g(ipjk2)),&
293  i_of(ijk2))
294  vgc2 = avg_x(avg_y_n(v_g(jm_of(ijk1)),v_g(ijk1)),&
295  avg_y_n(v_g(jm_of(ipjkp2)),v_g(ipjkp2)),&
296  i_of(ijk1))
297  vgc = avg_z(vgc1, vgc2, k_of(ijk2))
298  wgc = avg_x(w_g(ijk2), w_g(ipjk2),i_of(ijk2))
299  uscm = avg_z(u_s(ijk2, m), u_s(ijk1, m), k_of(ijk2))
300  vscm1= avg_x(avg_y_n(v_s(jm_of(ijk2),m),v_s(ijk2,m)),&
301  avg_y_n(v_s(jm_of(ipjk2),m),v_s(ipjk2,m)),&
302  i_of(ijk2))
303  vscm2= avg_x(avg_y_n(v_s(jm_of(ijk1),m),v_s(ijk1,m)),&
304  avg_y_n(v_s(jm_of(ipjkp2),m),v_s(ipjkp2,m)),&
305  i_of(ijk1))
306  vscm = avg_z(vscm1, vscm2, k_of(ijk2))
307  wscm = avg_x(w_s(ijk2, m), w_s(ipjk2, m),i_of(ijk2))
308  vels = uscm
309 
310  ELSE
311  WRITE(line,'(A, A)') 'Error: Unknown FCELL'
312  CALL write_error('CALC_GRBDRY', line, 1)
313  CALL exitmpi(mype)
314  ENDIF
315 
316 !-----------------------------------------------
317 ! Calculations for V momentum equation
318  ELSEIF (com .EQ. 'V')THEN
319 
320  ijk2n = north_of(ijk2)
321  ijpk2 = jp_of(ijk2)
322 
323  epg_avg = avg_y(ep_g(ijk2), ep_g(ijk2n), j_of(ijk2))
324  ep_star_avg = avg_y(ep_star_array(ijk2), ep_star_array(ijk2n), j_of(ijk2))
325  mu_g_avg = avg_y(mu_g(ijk2), mu_g(ijk2n), j_of(ijk2))
326  ro_g_avg = avg_y(ro_g(ijk2), ro_g(ijk2n), j_of(ijk2))
327 
328  k_12_avg = zero
329  tau_12_avg = zero
330  tau_1_avg = zero
331  IF(kt_type_enum == simonin_1996 .OR. &
332  kt_type_enum == ahmadi_1995) THEN
333  tau_1_avg = avg_y(tau_1(ijk2), tau_1(ijk2n), j_of(ijk2))
334 ! Simonin only:
335  k_12_avg = avg_y(k_12(ijk2), k_12(ijk2n), j_of(ijk2))
336  tau_12_avg = avg_y(tau_12(ijk2), tau_12(ijk2n), j_of(ijk2))
337  ENDIF
338 
339 
340  g0eps_avg = zero
341  DO mm = 1, smax
342  g0(mm) = g_0avg(ijk2, ijk2n, 'Y', j_of(ijk2), m, mm)
343  eps_avg(mm) = avg_y(ep_s(ijk2, mm), ep_s(ijk2n, mm), j_of(ijk2))
344  dp_avg(mm) = avg_y(d_p(ijk2,mm), d_p(ijk2n,mm), j_of(ijk2))
345  g0eps_avg = g0eps_avg + g_0avg(ijk2, ijk2n, 'Y', j_of(ijk2), m, mm) &
346  * avg_y(ep_s(ijk2, mm), ep_s(ijk2n, mm), j_of(ijk2))
347  ros_avg(mm) = avg_y(ro_s(ijk2, mm), ro_s(ijk2n, mm), j_of(ijk2))
348  th_avg(mm) = avg_y(theta_m(ijk2,mm), theta_m(ijk2n,mm), j_of(ijk2))
349  IF(th_avg(mm) < zero) th_avg(mm) = smalltheta
350  ENDDO
351 
352  wvels = bc_vw_s(l, m)
353 
354  IF(fcell .EQ. 'T')THEN
355  ijpkm2 = km_of(ijpk2)
356 
357  DO mm = 1, smax
358  avgx1 = avg_z(theta_m(ijk1,mm), theta_m(ijk2,mm), k_of(ijk1))
359  avgx2 = avg_z(theta_m(ijpkm2,mm), theta_m(ijpk2,mm), k_of(ijpkm2))
360  IF(avgx1 < zero .AND. avgx2 > zero) avgx1 = avgx2
361  IF(avgx2 < zero .AND. avgx1 > zero) avgx2 = avgx1
362  IF(avgx1 < zero .AND. avgx2 < zero) THEN
363  th_avg(mm) = smalltheta
364  ELSE
365  th_avg(mm) = avg_y(avgx1, avgx2, j_of(ijk1))
366  ENDIF
367  ENDDO
368 
369 ! Calculate velocity components at i,j+1/2,k+1/2 (relative to IJK1)
370  ugc1 = avg_x_e(&
371  avg_y(u_g(im_of(ijk1)), u_g(im_of(ijpkm2)),&
372  j_of(im_of(ijk1))),&
373  avg_y(u_g(ijk1), u_g(ijpkm2),j_of(ijk1)),&
374  i_of(ijk1))
375  ugc2 = avg_x_e(&
376  avg_y(u_g(im_of(ijk2)), u_g(im_of(ijpk2)),&
377  j_of(im_of(ijk2))),&
378  avg_y(u_g(ijk2), u_g(ijpk2),j_of(ijk2)),&
379  i_of(ijk2))
380  ugc = avg_z(ugc1, ugc2, k_of(ijk1))
381  vgc = avg_z(v_g(ijk1), v_g(ijk2),k_of(ijk1))
382  wgc = avg_y(w_g(ijk1), w_g(ijpkm2), j_of(ijk1))
383  uscm1= avg_x_e(&
384  avg_y(u_s(im_of(ijk1),m),u_s(im_of(ijpkm2),m),&
385  j_of(im_of(ijk1))),&
386  avg_y(u_s(ijk1,m), u_s(ijpkm2,m),j_of(ijk1)),&
387  i_of(ijk1))
388  uscm2= avg_x_e(&
389  avg_y(u_s(im_of(ijk2),m),u_s(im_of(ijpk2),m),&
390  j_of(im_of(ijk2))),&
391  avg_y(u_s(ijk2,m), u_s(ijpk2,m),j_of(ijk2)),&
392  i_of(ijk2))
393  uscm = avg_z(uscm1, uscm2, k_of(ijk1))
394  vscm = avg_z(v_s(ijk1,m), v_s(ijk2,m),k_of(ijk1))
395  wscm = avg_y(w_s(ijk1,m), w_s(ijpkm2,m), j_of(ijk1))
396  vels = vscm
397 
398  ELSEIF(fcell .EQ. 'B')THEN
399  ijpkp2 = kp_of(ijpk2)
400 
401  DO mm = 1, smax
402  avgx1 = avg_z(theta_m(ijk2,mm), theta_m(ijk1,mm), k_of(ijk2))
403  avgx2 = avg_z(theta_m(ijpk2,mm), theta_m(ijpkp2,mm), k_of(ijpk2))
404  IF(avgx1 < zero .AND. avgx2 > zero) avgx1 = avgx2
405  IF(avgx2 < zero .AND. avgx1 > zero) avgx2 = avgx1
406  IF(avgx1 < zero .AND. avgx2 < zero) THEN
407  th_avg(mm) = smalltheta
408  ELSE
409  th_avg(mm) = avg_y(avgx1, avgx2, j_of(ijk2))
410  ENDIF
411  ENDDO
412 
413 ! Calculate velocity components at i,j+1/2,k+1/2 (relative to IJK2)
414  ugc1 = avg_x_e(&
415  avg_y(u_g(im_of(ijk1)), u_g(im_of(ijpkp2)),&
416  j_of(im_of(ijk1))),&
417  avg_y(u_g(ijk1), u_g(ijpkp2),j_of(ijk1)),&
418  i_of(ijk1))
419  ugc2 = avg_x_e(&
420  avg_y(u_g(im_of(ijk2)), u_g(im_of(ijpk2)),&
421  j_of(im_of(ijk2))),&
422  avg_y(u_g(ijk2), u_g(ijpk2),j_of(ijk2)),&
423  i_of(ijk2))
424  ugc = avg_z(ugc2, ugc1, k_of(ijk2))
425  vgc = avg_z(v_g(ijk2), v_g(ijk1),k_of(ijk2))
426  wgc = avg_y(w_g(ijk2), w_g(ijpk2), j_of(ijk2))
427  uscm1= avg_x_e(&
428  avg_y(u_s(im_of(ijk1),m),u_s(im_of(ijpkp2),m),&
429  j_of(im_of(ijk1))),&
430  avg_y(u_s(ijk1,m), u_s(ijpkp2,m),j_of(ijk1)),&
431  i_of(ijk1))
432  uscm2= avg_x_e(&
433  avg_y(u_s(im_of(ijk2),m),u_s(im_of(ijpk2),m),&
434  j_of(im_of(ijk2))),&
435  avg_y(u_s(ijk2,m), u_s(ijpk2,m),j_of(ijk2)),&
436  i_of(ijk2))
437  uscm = avg_z(uscm2, uscm1, k_of(ijk2))
438  vscm = avg_z(v_s(ijk2,m), v_s(ijk1,m),k_of(ijk2))
439  wscm = avg_y(w_s(ijk2,m), w_s(ijpk2,m), j_of(ijk2))
440  vels = vscm
441 
442  ELSEIF(fcell .EQ. 'E')THEN
443  imjpk2= im_of(ijpk2)
444 
445  DO mm = 1, smax
446  avgx1 = avg_x(theta_m(ijk1,mm),theta_m(ijk2,mm),i_of(ijk1))
447  avgx2 = avg_x(theta_m(imjpk2,mm),theta_m(ijpk2,mm),i_of(imjpk2))
448  IF(avgx1 < zero .AND. avgx2 > zero) avgx1 = avgx2
449  IF(avgx2 < zero .AND. avgx1 > zero) avgx2 = avgx1
450  IF(avgx1 < zero .AND. avgx2 < zero) THEN
451  th_avg(mm) = smalltheta
452  ELSE
453  th_avg(mm) = avg_y(avgx1, avgx2, j_of(ijk1))
454  ENDIF
455  ENDDO
456 
457 ! Calculate velocity components at i+1/2,j+1/2,k relative to IJK1
458  ugc = avg_y(u_g(ijk1), u_g(imjpk2), j_of(ijk1))
459  vgc = avg_x(v_g(ijk1), v_g(ijk2), i_of(ijk1))
460  wgc1 = avg_y(avg_z_t(w_g(km_of(ijk1)), w_g(ijk1)),&
461  avg_z_t(w_g(km_of(imjpk2)), w_g(imjpk2)),&
462  j_of(ijk1))
463  wgc2 = avg_y(avg_z_t(w_g(km_of(ijk2)), w_g(ijk2)),&
464  avg_z_t(w_g(km_of(ijpk2)), w_g(ijpk2)),&
465  j_of(ijk2))
466  wgc = avg_x(wgc1, wgc2, i_of(ijk1))
467  uscm = avg_y(u_s(ijk1,m), u_s(imjpk2,m), j_of(ijk1))
468  vscm = avg_x(v_s(ijk1, m), v_s(ijk2, m), i_of(ijk1))
469  wscm1= avg_y(avg_z_t(w_s(km_of(ijk1),m), w_s(ijk1,m)),&
470  avg_z_t(w_s(km_of(imjpk2),m), w_s(imjpk2,m)),&
471  j_of(ijk1))
472  wscm2 = avg_y(avg_z_t(w_s(km_of(ijk2),m), w_s(ijk2,m)),&
473  avg_z_t(w_s(km_of(ijpk2),m), w_s(ijpk2,m)),&
474  j_of(ijk2))
475  wscm = avg_x(wscm1, wscm2, i_of(ijk1))
476  vels = vscm
477 
478  ELSEIF(fcell .EQ. 'W')THEN
479  ipjpk2= ip_of(ijpk2)
480 
481  DO mm = 1, smax
482  avgx1 = avg_x(theta_m(ijk2,mm),theta_m(ijk1,mm),i_of(ijk2))
483  avgx2 = avg_x(theta_m(ijpk2,mm),theta_m(ipjpk2,mm),i_of(ijpk2))
484  IF(avgx1 < zero .AND. avgx2 > zero) avgx1 = avgx2
485  IF(avgx2 < zero .AND. avgx1 > zero) avgx2 = avgx1
486  IF(avgx1 < zero .AND. avgx2 < zero) THEN
487  th_avg(mm) = smalltheta
488  ELSE
489  th_avg(mm) = avg_y(avgx1, avgx2, j_of(ijk2))
490  ENDIF
491  ENDDO
492 
493 ! Calculate velocity components at i+1/2,j+1/2,k relative to IJK2
494  ugc = avg_y(u_g(ijk2), u_g(ijpk2), j_of(ijk2))
495  vgc = avg_x(v_g(ijk2), v_g(ijk1), i_of(ijk2))
496  wgc1 = avg_y(avg_z_t(w_g(km_of(ijk1)), w_g(ijk1)),&
497  avg_z_t(w_g(km_of(ipjpk2)), w_g(ipjpk2)),&
498  j_of(ijk1))
499  wgc2 = avg_y(avg_z_t(w_g(km_of(ijk2)), w_g(ijk2)),&
500  avg_z_t(w_g(km_of(ijpk2)), w_g(ijpk2)),&
501  j_of(ijk2))
502  wgc = avg_x(wgc2, wgc1, i_of(ijk2))
503  uscm = avg_y(u_s(ijk2,m), u_s(ijpk2,m), j_of(ijk2))
504  vscm = avg_x(v_s(ijk2, m), v_s(ijk1, m), i_of(ijk2))
505  wscm1= avg_y(avg_z_t(w_s(km_of(ijk1),m), w_s(ijk1,m)),&
506  avg_z_t(w_s(km_of(ipjpk2),m), w_s(ipjpk2,m)),&
507  j_of(ijk1))
508  wscm2 = avg_y(avg_z_t(w_s(km_of(ijk2),m), w_s(ijk2,m)),&
509  avg_z_t(w_s(km_of(ijpk2),m), w_s(ijpk2,m)),&
510  j_of(ijk2))
511  wscm = avg_x(wscm2, wscm1, i_of(ijk2))
512  vels = vscm
513 
514  ELSE
515  WRITE(line,'(A, A)') 'Error: Unknown FCELL'
516  CALL write_error('CALC_GRBDRY', line, 1)
517  CALL exitmpi(mype)
518  ENDIF
519 
520 
521 !-----------------------------------------------
522 ! Calculations for W momentum equation
523  ELSEIF (com .EQ. 'W')THEN
524  ijk2t = top_of(ijk2)
525  ijkp2 = kp_of(ijk2)
526 
527  epg_avg = avg_z(ep_g(ijk2), ep_g(ijk2t), k_of(ijk2))
528  mu_g_avg = avg_z(mu_g(ijk2), mu_g(ijk2t), k_of(ijk2))
529  ro_g_avg = avg_z(ro_g(ijk2), ro_g(ijk2t), k_of(ijk2))
530  ep_star_avg = avg_z(ep_star_array(ijk2), ep_star_array(ijk2t), k_of(ijk2))
531 
532  k_12_avg = zero
533  tau_12_avg = zero
534  tau_1_avg = zero
535  IF(kt_type_enum == simonin_1996 .OR. &
536  kt_type_enum == ahmadi_1995) THEN
537  tau_1_avg = avg_z(tau_1(ijk2), tau_1(ijk2t), k_of(ijk2))
538 ! Simonin only:
539  k_12_avg = avg_z(k_12(ijk2), k_12(ijk2t), k_of(ijk2))
540  tau_12_avg = avg_z(tau_12(ijk2), tau_12(ijk2t), k_of(ijk2))
541  ENDIF
542 
543  g0eps_avg = zero
544  DO mm = 1, smax
545  g0(mm) = g_0avg(ijk2, ijk2t, 'Z', k_of(ijk2), m, mm)
546  eps_avg(mm) = avg_z(ep_s(ijk2,mm), ep_s(ijk2t,mm), k_of(ijk2))
547  dp_avg(mm) = avg_z(d_p(ijk2,mm), d_p(ijk2t,mm), k_of(ijk2))
548  g0eps_avg = g0eps_avg + g_0avg(ijk2, ijk2t, 'Z', k_of(ijk2), m, mm) &
549  * avg_z(ep_s(ijk2, mm), ep_s(ijk2t, mm), k_of(ijk2))
550  ros_avg(mm) = avg_z(ro_s(ijk2,mm), ro_s(ijk2t,mm), k_of(ijk2))
551  th_avg(mm) = avg_z(theta_m(ijk2,mm), theta_m(ijk2t,mm), k_of(ijk2))
552  IF(th_avg(mm) < zero) th_avg(mm) = smalltheta
553  ENDDO
554 
555  wvels = bc_ww_s(l,m)
556 
557  IF(fcell .EQ. 'N')THEN
558  ijmkp2 = jm_of(ijkp2)
559 
560  DO mm = 1, smax
561  avgx1 = avg_z(theta_m(ijk1,mm), theta_m(ijmkp2,mm), k_of(ijk1))
562  avgx2 = avg_z(theta_m(ijk2,mm), theta_m(ijkp2,mm), k_of(ijk2))
563  IF(avgx1 < zero .AND. avgx2 > zero) avgx1 = avgx2
564  IF(avgx2 < zero .AND. avgx1 > zero) avgx2 = avgx1
565  IF(avgx1 < zero .AND. avgx2 < zero) THEN
566  th_avg(mm) = smalltheta
567  ELSE
568  th_avg(mm) = avg_y(avgx1, avgx2, j_of(ijk1))
569  ENDIF
570  ENDDO
571 
572 ! Calculate velocity components at i,j+1/2,k+1/2 (relative to IJK1)
573  ugc1 = avg_x_e(&
574  avg_z(u_g(im_of(ijk1)), u_g(im_of(ijmkp2)),&
575  k_of(im_of(ijk1)) ),&
576  avg_z(u_g(ijk1), u_g(ijmkp2), k_of(ijk1)),&
577  i_of(ijk1))
578  ugc2 = avg_x_e(&
579  avg_z(u_g(im_of(ijk2)), u_g(im_of(ijkp2)),&
580  k_of(im_of(ijk2))),&
581  avg_z(u_g(ijk2), u_g(ijkp2), k_of(ijk2)),&
582  i_of(ijk2))
583  ugc = avg_y(ugc1, ugc2, j_of(ijk1))
584  vgc = avg_z(v_g(ijk1), v_g(ijmkp2),k_of(ijk1))
585  wgc = avg_y(w_g(ijk1), w_g(ijk2), j_of(ijk1))
586  uscm1= avg_x_e(&
587  avg_z(u_s(im_of(ijk1),m),u_s(im_of(ijmkp2),m),&
588  k_of(im_of(ijk1))),&
589  avg_z(u_s(ijk1,m), u_s(ijmkp2,m),k_of(ijk1)),&
590  i_of(ijk1))
591  uscm2= avg_x_e(&
592  avg_z(u_s(im_of(ijk2),m),u_s(im_of(ijkp2),m),&
593  k_of(im_of(ijk2))),&
594  avg_z(u_s(ijk2,m), u_s(ijkp2,m),k_of(ijk2)),&
595  i_of(ijk2))
596  uscm = avg_y(uscm1, uscm2, j_of(ijk1))
597  vscm = avg_z(v_s(ijk1,m), v_s(ijmkp2,m),k_of(ijk1))
598  wscm = avg_y(w_s(ijk1,m), w_s(ijk2,m), j_of(ijk1))
599  vels = wscm
600 
601  ELSEIF(fcell .EQ. 'S')THEN
602  ijpkp2 = jp_of(ijkp2)
603 
604  DO mm = 1, smax
605  avgx1 = avg_z(theta_m(ijk2,mm), theta_m(ijkp2,mm), k_of(ijk2))
606  avgx2 = avg_z(theta_m(ijk1,mm), theta_m(ijpkp2,mm), k_of(ijk1))
607  IF(avgx1 < zero .AND. avgx2 > zero) avgx1 = avgx2
608  IF(avgx2 < zero .AND. avgx1 > zero) avgx2 = avgx1
609  IF(avgx1 < zero .AND. avgx2 < zero) THEN
610  th_avg(mm) = smalltheta
611  ELSE
612  th_avg(mm) = avg_y(avgx1, avgx2, j_of(ijk2))
613  ENDIF
614  ENDDO
615 
616 ! Calculate velocity components at i,j+1/2,k+1/2 (relative to IJK2)
617  ugc1 = avg_x_e(&
618  avg_z(u_g(im_of(ijk1)), u_g(im_of(ijpkp2)),&
619  k_of(im_of(ijk1))),&
620  avg_z(u_g(ijk1), u_g(ijpkp2), k_of(ijk1)),&
621  i_of(ijk1))
622  ugc2 = avg_x_e(&
623  avg_z(u_g(im_of(ijk2)), u_g(im_of(ijkp2)),&
624  k_of(im_of(ijk2))),&
625  avg_z(u_g(ijk2), u_g(ijkp2), k_of(ijk2)),&
626  i_of(ijk2))
627  ugc = avg_y(ugc2, ugc1, j_of(ijk2))
628  vgc = avg_z(v_g(ijk2), v_g(ijkp2),k_of(ijk2))
629  wgc = avg_y(w_g(ijk2), w_g(ijk1), j_of(ijk2))
630  uscm1= avg_x_e(&
631  avg_z(u_s(im_of(ijk1),m),u_s(im_of(ijpkp2),m),&
632  k_of(im_of(ijk1))),&
633  avg_z(u_s(ijk1,m), u_s(ijpkp2,m),k_of(ijk1)),&
634  i_of(ijk1))
635  uscm2= avg_x_e(&
636  avg_z(u_s(im_of(ijk2),m),u_s(im_of(ijkp2),m),&
637  k_of(im_of(ijk2))),&
638  avg_z(u_s(ijk2,m), u_s(ijkp2,m),k_of(ijk2)),&
639  i_of(ijk2))
640  uscm = avg_y(uscm2, uscm1, j_of(ijk2))
641  vscm = avg_z(v_s(ijk2,m), v_s(ijkp2,m),k_of(ijk2))
642  wscm = avg_y(w_s(ijk2,m), w_s(ijk1,m), j_of(ijk2))
643  vels = wscm
644 
645  ELSEIF(fcell .EQ. 'E')THEN
646  imjkp2 = im_of(ijkp2)
647 
648  DO mm = 1, smax
649  avgx1 = avg_x(theta_m(ijk1,mm),theta_m(ijk2,mm),i_of(ijk1))
650  avgx2 = avg_x(theta_m(imjkp2,mm),theta_m(ijkp2,mm),i_of(imjkp2))
651  IF(avgx1 < zero .AND. avgx2 > zero) avgx1 = avgx2
652  IF(avgx2 < zero .AND. avgx1 > zero) avgx2 = avgx1
653  IF(avgx1 < zero .AND. avgx2 < zero) THEN
654  th_avg(mm) = smalltheta
655  ELSE
656  th_avg(mm) = avg_z(avgx1, avgx2, k_of(ijk1))
657  ENDIF
658  ENDDO
659 
660 ! Calculate velocity components at i+1/2,j,k+1/2 relative to IJK1
661  ugc = avg_z(u_g(ijk1), u_g(imjkp2), k_of(ijk1))
662  vgc1 = avg_z(avg_y_n(v_g(jm_of(ijk1)),v_g(ijk1)),&
663  avg_y_n(v_g(jm_of(imjkp2)),v_g(imjkp2)),&
664  k_of(ijk1))
665  vgc2 = avg_z(avg_y_n(v_g(jm_of(ijk2)),v_g(ijk2)),&
666  avg_y_n(v_g(jm_of(ijkp2)),v_g(ijkp2)),&
667  k_of(ijk2))
668  vgc = avg_x(vgc1,vgc2,i_of(ijk1))
669  wgc = avg_x(w_g(ijk1), w_g(ijk2),i_of(ijk1))
670  uscm = avg_z(u_s(ijk1,m), u_s(imjkp2,m), k_of(ijk1))
671  vscm1= avg_z(avg_y_n(v_s(jm_of(ijk1),m),v_s(ijk1,m)),&
672  avg_y_n(v_s(jm_of(imjkp2),m),v_s(imjkp2,m)),&
673  k_of(ijk1))
674  vscm2= avg_z(avg_y_n(v_s(jm_of(ijk2),m),v_s(ijk2,m)),&
675  avg_y_n(v_s(jm_of(ijkp2),m),v_s(ijkp2,m)),&
676  k_of(ijk2))
677  vscm = avg_x(vscm1,vscm2,i_of(ijk1))
678  wscm = avg_x(w_s(ijk1,m), w_s(ijk2,m), i_of(ijk1))
679  vels = wscm
680 
681  ELSEIF(fcell .EQ. 'W')THEN
682  ipjkp2= ip_of(ijkp2)
683 
684  DO mm = 1, smax
685  avgx1 = avg_x(theta_m(ijk2,mm),theta_m(ijk1,mm),i_of(ijk2))
686  avgx2 = avg_x(theta_m(ijkp2,mm),theta_m(ipjkp2,mm),i_of(ijkp2))
687  IF(avgx1 < zero .AND. avgx2 > zero) avgx1 = avgx2
688  IF(avgx2 < zero .AND. avgx1 > zero) avgx2 = avgx1
689  IF(avgx1 < zero .AND. avgx2 < zero) THEN
690  th_avg(mm) = smalltheta
691  ELSE
692  th_avg(mm) = avg_z(avgx1, avgx2, k_of(ijk2))
693  ENDIF
694  ENDDO
695 
696 ! Calculate velocity components at i+1/2,j,k+1/2 relative to IJK2
697  ugc = avg_z(u_g(ijk2), u_g(ijkp2), k_of(ijk2))
698  vgc1 = avg_z(avg_y_n(v_g(jm_of(ijk1)),v_g(ijk1)),&
699  avg_y_n(v_g(jm_of(ipjkp2)),v_g(ipjkp2)),&
700  k_of(ijk1))
701  vgc2 = avg_z(avg_y_n(v_g(jm_of(ijk2)),v_g(ijk2)),&
702  avg_y_n(v_g(jm_of(ijkp2)),v_g(ijkp2)),&
703  k_of(ijk2))
704  vgc = avg_x(vgc2,vgc1,i_of(ijk2))
705  wgc = avg_x(w_g(ijk2), w_g(ijk1),i_of(ijk2))
706  uscm = avg_z(u_s(ijk2,m), u_s(ijkp2,m), k_of(ijk2))
707  vscm1= avg_z(avg_y_n(v_s(jm_of(ijk1),m),v_s(ijk1,m)),&
708  avg_y_n(v_s(jm_of(ipjkp2),m),v_s(ipjkp2,m)),&
709  k_of(ijk1))
710  vscm2= avg_z(avg_y_n(v_s(jm_of(ijk2),m),v_s(ijk2,m)),&
711  avg_y_n(v_s(jm_of(ijkp2),m),v_s(ijkp2,m)),&
712  k_of(ijk2))
713  vscm = avg_x(vscm2,vscm1,i_of(ijk2))
714  wscm = avg_x(w_s(ijk2,m), w_s(ijk1,m), i_of(ijk2))
715  vels = wscm
716 
717  ELSE
718  WRITE(line,'(A, A)') 'Error: Unknown FCELL'
719  CALL write_error('CALC_GRBDRY', line, 1)
720  CALL exitmpi(mype)
721  ENDIF
722 
723 
724  ELSE
725  WRITE(line,'(A, A)') 'Error: Unknown COM'
726  CALL write_error('CALC_GRBDRY', line, 1)
727  CALL exitmpi(mype)
728  ENDIF
729 
730 ! magnitude of gas-solids relative velocity
731  vrel =&
732  dsqrt( (ugc - uscm)**2 + (vgc - vscm)**2 + (wgc - wscm)**2 )
733 
734 ! slip velocity for use in Jenkins bc (sof)
735  vslip= dsqrt( (uscm-bc_uw_s(l,m))**2 + (vscm-bc_vw_s(l,m))**2 &
736  + (wscm-bc_ww_s(l,m))**2 )
737 
738  IF (friction .AND. (one-ep_g(ijk2))>eps_f_min) THEN
739  CALL calc_s_ddot_s(ijk1, ijk2, fcell, com, m, del_dot_u,&
740  s_ddot_s, s_dd)
741 
742  CALL calc_gw_hw_cw(g0(m), eps_avg(m), epg_avg, ep_star_avg, &
743  g0eps_avg, th_avg(m), mu_g_avg, ro_g_avg, ros_avg(m), &
744  dp_avg(m), k_12_avg, tau_12_avg, tau_1_avg, vrel, vslip,&
745  del_dot_u, s_ddot_s, s_dd, vels, wvels, m, gw, hw, cw)
746  ELSE
747  gw = 1d0
748  hw = f_hw(g0, eps_avg, epg_avg, ep_star_avg, &
749  g0eps_avg, th_avg, mu_g_avg, ro_g_avg, ros_avg, &
750  dp_avg, k_12_avg, tau_12_avg, tau_1_avg, vrel, vslip, m)
751  cw = hw*wvels
752  ENDIF
753 
754 
755  RETURN
756  END SUBROUTINE calc_grbdry
757 
758 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
759 ! C
760 ! Function: F_HW C
761 ! Purpose: Function for hw C
762 ! C
763 ! Literature/Document References: C
764 ! See calc_mu_s.f for references on kinetic theory models C
765 ! Johnson, P. C., and Jackson, R., Frictional-collisional C
766 ! constitutive relations for granular materials, with C
767 ! application to plane shearing, Journal of Fluid Mechanics, C
768 ! 1987, 176, pp. 67-93 C
769 ! Jenkins, J. T., and Louge, M. Y., On the flux of fluctuating C
770 ! energy in a collisional grain flow at a flat frictional wall, C
771 ! Physics of Fluids, 1997, 9(10), pp. 2835-2840 C
772 ! C
773 ! C
774 ! Additional Notes: C
775 ! The JJ granular momentum BC is written as: C
776 ! usw/|usw|.tau.n+C+F=0 C
777 ! tau = the stress tensor (collisional+frictional) C
778 ! n = outward normal vector C
779 ! usw = slip velocity with wall C
780 ! C = collisional force C
781 ! F = frictional force C
782 ! C
783 ! The granular momentum BC is written as the normal vector dot the C
784 ! stress tensor. Besides the gradient in velocity of phase M, the C
785 ! stress tensor expression may contain several additional terms C
786 ! that would need to be accounted for when satisfying the BC. These C
787 ! modifications have NOT been rigorously addressed. C
788 ! C
789 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
790  DOUBLE PRECISION FUNCTION f_hw(g0,EPS,EPG, ep_star_avg, &
791  g0eps_avg,th,mu_g_avg,ro_g_avg,ros_avg, &
792  dp_avg,k_12_avg, tau_12_avg, tau_1_avg, &
793  vrel, vslip, m)
794 
795 !-----------------------------------------------
796 ! Modules
797 !-----------------------------------------------
798  USE param
799  USE param1
800  USE constant
801  USE kintheory, only: epm, k_phi, s_star
802  USE physprop
803  USE run
804  USE fldvar
805  USE mpi_utility
806  IMPLICIT NONE
807 !-----------------------------------------------
808 ! Dummy Arguments
809 !-----------------------------------------------
810 ! Radial distribution function of solids phase M with each
811 ! other solids phase
812  DOUBLE PRECISION, INTENT(IN) :: g0(dimension_m)
813 ! Average solids volume fraction of each solids phase
814  DOUBLE PRECISION, INTENT(IN) :: EPS(dimension_m)
815 ! Average solids and gas volume fraction
816  DOUBLE PRECISION, INTENT(IN) :: EPG, ep_star_avg
817 ! Sum of eps*G_0
818  DOUBLE PRECISION, INTENT(INOUT) :: g0EPs_avg
819 ! Average theta_m
820  DOUBLE PRECISION, INTENT(INOUT) :: TH (dimension_m)
821 ! Average gas viscosity
822  DOUBLE PRECISION, INTENT(IN) :: Mu_g_avg
823 ! Average gas density
824  DOUBLE PRECISION, INTENT(IN) :: RO_g_avg
825 ! Average solids density
826  DOUBLE PRECISION, INTENT(IN) :: ROS_avg(dimension_m)
827 ! Average particle diameter of each solids phase
828  DOUBLE PRECISION, INTENT(IN) :: DP_avg(dimension_m)
829 ! Average cross-correlation K_12 and lagrangian integral time-scale
830  DOUBLE PRECISION, INTENT(IN) :: K_12_avg, Tau_12_avg, Tau_1_avg
831 ! Magnitude of slip velocity between two phases
832  DOUBLE PRECISION, INTENT(IN) :: VREL
833 ! Slip velocity between wall and particles
834  DOUBLE PRECISION, INTENT(IN) :: VSLIP
835 ! Solids phase index
836  INTEGER, INTENT(IN) :: M
837 !-----------------------------------------------
838 ! Local Variables
839 !-----------------------------------------------
840 ! Solids phase index
841  INTEGER :: LL
842 ! Coefficient of 2nd term
843  DOUBLE PRECISION :: F_2
844 ! Coefficient of 1st term
845  DOUBLE PRECISION :: Mu_s
846 ! Viscosity
847  DOUBLE PRECISION :: Mu
848 ! Bulk viscosity
849  DOUBLE PRECISION :: Mu_b
850 ! Viscosity corrected for interstitial fluid effects
851  DOUBLE PRECISION :: Mu_star
852 ! Solids pressure
853  DOUBLE PRECISION :: Ps
854 ! Reynolds number based on slip velocity
855  DOUBLE PRECISION :: Re_g
856 ! Friction Factor in drag coefficient
857  DOUBLE PRECISION :: C_d
858 ! Drag Coefficient
859  DOUBLE PRECISION :: Beta, DgA
860 ! Constants in Simonin or Ahmadi model
861  DOUBLE PRECISION :: Sigma_c, Tau_2_c, Tau_12_st, Nu_t
862  DOUBLE PRECISION :: Tau_2, zeta_c_2, MU_2_T_Kin, Mu_2_Col
863  DOUBLE PRECISION :: Tmp_Ahmadi_Const
864 ! Variables for Iddir & Arastoopour model
865  DOUBLE PRECISION :: MU_sM_sum, MU_s_MM, MU_s_LM, MU_sM_ip, MU_common_term,&
866  MU_sM_LM
867  DOUBLE PRECISION :: M_PM, M_PL, MPSUM, NU_PL, NU_PM, D_PM, D_PL, DPSUMo2
868  DOUBLE PRECISION :: Ap_lm, Dp_lm, R1p_lm, Bp_lm
869 ! Variables for Garzo and Dufty model
870  DOUBLE PRECISION :: c_star, zeta0_star, nu_eta_star, press_star, &
871  gamma_star, eta_k_star, eta_star, eta0
872 ! Variables for GTSH theory
873  DOUBLE PRECISION :: Chi, RdissP, Re_T, Rdiss, GamaAvg, A2_gtshW, &
874  zeta_star, nu0, nuN, NuK, EDT_s, zeta_avg, etaK, &
875  Mu_b_avg, mu2_0, mu4_0, mu4_1
876 !-----------------------------------------------
877 ! Functions
878 !-----------------------------------------------
879 ! Variable specularity coefficient
880  DOUBLE PRECISION :: PHIP_JJ
881 !-----------------------------------------------
882 
883 ! This is done here similar to bc_theta to avoid small negative values of
884 ! Theta coming most probably from linear solver
885  IF(th(m) .LE. zero)THEN
886  th(m) = 1d-8
887  IF (mype.eq.pe_io) THEN
888  WRITE(*,*) &
889  'Warning: Negative granular temp at wall set to 1e-8'
890 ! CALL WRITE_ERROR('THETA_HW_CW', LINE, 1)
891  ENDIF
892  ENDIF
893 
894 
895 ! common variables
896  d_pm = dp_avg(m)
897  m_pm = (pi/6.d0)*(d_pm**3)*ros_avg(m)
898  nu_pm = (eps(m)*ros_avg(m))/m_pm
899 
900 ! This is from Wen-Yu correlation, you can put here your own single particle drag
901  re_g = epg*ro_g_avg*d_pm*vrel/mu_g_avg
902  IF (re_g .lt. 1000.d0) THEN
903  c_d = (24.d0/(re_g+small_number))*(one + 0.15d0 * re_g**0.687d0)
904  ELSE
905  c_d = 0.44d0
906  ENDIF
907  dga = 0.75d0*c_d*ro_g_avg*epg*vrel/(d_pm*epg**(2.65d0))
908  IF(vrel == zero) dga = large_number
909  beta = eps(m)*dga !this is equivalent to F_gs(ijk,m)
910 
911 
912 ! Defining viscosity according to each KT for later use in JJ BC
913 ! Also defining solids pressure according to each KT for use if JENKINS
914 ! -------------------------------------------------------------------
915  SELECT CASE (kt_type_enum)
916  CASE (lun_1984)
917  mu = (5d0*dsqrt(pi*th(m))*dp_avg(m)*ros_avg(m))/96d0
918  mu_b = (256d0*mu*eps(m)*g0eps_avg)/(5d0*pi)
919 
920  IF(switch == zero .OR. ro_g_avg == zero)THEN
921  mu_star = mu
922  ELSEIF(th(m) .LT. small_number)THEN
923  mu_star = zero
924  ELSE
925  mu_star = ros_avg(m)*eps(m)* g0(m)*th(m)* mu/ &
926  (ros_avg(m)*g0eps_avg*th(m) + 2.0d0*dga/ros_avg(m)* mu)
927  ENDIF
928 
929  mu_s = ((2d0+alpha)/3d0)*((mu_star/(eta*(2d0-eta)*&
930  g0(m)))*(one+1.6d0*eta*g0eps_avg)*&
931  (one+1.6d0*eta*(3d0*eta-2d0)*&
932  g0eps_avg)+(0.6d0*mu_b*eta))
933 
934 ! defining granular pressure (for Jenkins BC)
935  ps = ros_avg(m)*eps(m)*th(m)*(one+4.d0*eta*g0eps_avg)
936 
937 
938  CASE (simonin_1996)
939 ! particle relaxation time
940  tau_12_st = ros_avg(m)/(dga+small_number)
941 
942  sigma_c = (one+ c_e)*(3.d0-c_e)/5.d0
943  tau_2_c = dp_avg(m)/(6.d0*eps(m)*g0(m)*dsqrt(16.d0*(th(m)+small_number)/pi))
944  zeta_c_2= 2.d0/5.d0*(one+ c_e)*(3.d0*c_e-one)
945  nu_t = tau_12_avg/tau_12_st
946  tau_2 = one/(2.d0/tau_12_st+sigma_c/tau_2_c)
947  mu_2_t_kin = (2.0d0/3.0d0*k_12_avg*nu_t + th(m) * &
948  (one+ zeta_c_2*eps(m)*g0(m)))*tau_2
949  mu_2_col = 8.d0/5.d0*eps(m)*g0(m)*eta* (mu_2_t_kin+ &
950  dp_avg(m)*dsqrt(th(m)/pi))
951  mu_s = eps(m)*ros_avg(m)*(mu_2_t_kin + mu_2_col)
952 
953 ! defining granular pressure (for Jenkins BC)
954  ps = ros_avg(m)*eps(m)*th(m)*(one+4.d0*eta*g0eps_avg)
955 
956 
957  CASE (ahmadi_1995)
958 ! particle relaxation time
959  tau_12_st = ros_avg(m)/(dga+small_number)
960 
961  IF(eps(m) < (one-ep_star_avg)) THEN
962  tmp_ahmadi_const = one/&
963  (one+ tau_1_avg/tau_12_st * (one-eps(m)/(one-ep_star_avg))**3)
964  ELSE
965  tmp_ahmadi_const = one
966  ENDIF
967  mu_s = tmp_ahmadi_const &
968  *0.1045d0*(one/g0(m)+3.2d0*eps(m)+12.1824d0*g0(m)*eps(m)*eps(m)) &
969  *dp_avg(m)*ros_avg(m)* dsqrt(th(m))
970 
971 ! defining granular pressure (for Jenkins BC)
972  ps = ros_avg(m)*eps(m)*th(m)*&
973  ((one + 4.0d0*g0eps_avg)+half*(one -c_e*c_e))
974 
975 
976  CASE(ia_2005)
977 ! Use original IA theory if SWITCH_IA is false
978  IF(.NOT. switch_ia) g0eps_avg = eps(m)*ros_avg(m)
979 
980  mu = (5.d0/96.d0)*d_pm* ros_avg(m)*dsqrt(pi*th(m)/m_pm)
981 
982  IF(switch == zero .OR. ro_g_avg == zero)THEN
983  mu_star = mu
984  ELSEIF(th(m) .LT. small_number)THEN
985  mu_star = zero
986  ELSE
987  mu_star = mu*eps(m)*g0(m)/ (g0eps_avg+ 2.0d0*dga*mu / &
988  (ros_avg(m)**2 *(th(m)/m_pm)))
989  ENDIF
990 
991  mu_s_mm = (mu_star/g0(m))*(1.d0+(4.d0/5.d0)*(1.d0+c_e)*g0eps_avg)**2
992  mu_sm_sum = zero
993 
994  DO ll = 1, smax
995  d_pl = dp_avg(ll)
996  m_pl = (pi/6.d0)*(d_pl**3.)*ros_avg(ll)
997  mpsum = m_pm + m_pl
998  dpsumo2 = (d_pm+d_pl)/2.d0
999  nu_pl = (eps(ll)*ros_avg(ll))/m_pl
1000 
1001  IF ( ll .eq. m) THEN
1002  ap_lm = mpsum/(2.d0)
1003  dp_lm = m_pl*m_pm/(2.d0*mpsum)
1004  r1p_lm = one/( ap_lm**1.5 * dp_lm**3 )
1005  mu_s_lm = dsqrt(pi)*(dpsumo2**4 / (48.d0*5.d0))*&
1006  g0(ll)*(m_pl*m_pm/mpsum)*(m_pl*m_pm/&
1007  mpsum)*((m_pl*m_pm)**1.5)*nu_pm*nu_pl*&
1008  (1.d0+c_e)*r1p_lm*dsqrt(th(m))
1009 
1010 ! solids phase 'viscosity' associated with the divergence
1011 ! of solids phase M
1012  mu_sm_ip = (mu_s_mm + mu_s_lm)
1013 
1014  ELSE
1015  ap_lm = (m_pm*th(ll)+m_pl*th(m))/&
1016  (2.d0)
1017  bp_lm = (m_pm*m_pl*(th(ll)-th(m) ))/&
1018  (2.d0*mpsum)
1019  dp_lm = (m_pl*m_pm*(m_pm*th(m)+m_pl*th(ll) ))/&
1020  (2.d0*mpsum*mpsum)
1021  r1p_lm = ( one/( ap_lm**1.5 * dp_lm**3 ) ) + &
1022  ( (9.d0*bp_lm*bp_lm)/( ap_lm**2.5 * dp_lm**4 ) )+&
1023  ( (30.d0*bp_lm**4)/( 2.d0*ap_lm**3.5 * dp_lm**5 ) )
1024  mu_common_term = dsqrt(pi)*(dpsumo2**4 / (48.d0*5.d0))*&
1025  g0(ll)*(m_pl*m_pm/mpsum)*(m_pl*m_pm/&
1026  mpsum)*((m_pl*m_pm)**1.5)*nu_pm*nu_pl*&
1027  (1.d0+c_e)*r1p_lm
1028  mu_sm_lm = mu_common_term*(th(m)*th(m)*th(ll)*th(ll)*th(ll) )
1029 
1030 ! solids phase 'viscosity' associated with the divergence
1031 ! of solids phase M
1032  mu_sm_ip = mu_sm_lm
1033 
1034  ENDIF
1035  mu_sm_sum = mu_sm_sum + mu_sm_ip
1036 
1037  ENDDO
1038 
1039 ! Find the term proportional to the gradient in velocity
1040 ! of phase M (viscosity in the Mth solids phase)
1041  mu_s = mu_sm_sum
1042 
1043 
1044  CASE(gd_1999)
1045 ! Pressure/Viscosity/Bulk Viscosity
1046 ! Note: k_boltz = M_PM
1047 
1048 ! Find pressure in the Mth solids phase
1049  press_star = 1.d0 + 2.d0*(1.d0+c_e)*eps(m)*g0(m)
1050 
1051 ! find bulk and shear viscosity
1052  c_star = 32.0d0*(1.0d0 - c_e)*(1.d0 - 2.0d0*c_e*c_e) &
1053  / (81.d0 - 17.d0*c_e + 30.d0*c_e*c_e*(1.0d0-c_e))
1054 
1055  zeta0_star = (5.d0/12.d0)*g0(m)*(1.d0 - c_e*c_e) &
1056  * (1.d0 + (3.d0/32.d0)*c_star)
1057 
1058  nu_eta_star = g0(m)*(1.d0 - 0.25d0*(1.d0-c_e)*(1.d0-c_e)) &
1059  * (1.d0-(c_star/64.d0))
1060 
1061  gamma_star = (4.d0/5.d0)*(32.d0/pi)*eps(m)*eps(m) &
1062  * g0(m)*(1.d0+c_e) * (1.d0 - (c_star/32.d0))
1063 
1064  eta_k_star = (1.d0 - (2.d0/5.d0)*(1.d0+c_e)*(1.d0-3.d0*c_e) &
1065  * eps(m)*g0(m) ) / (nu_eta_star - 0.5d0*zeta0_star)
1066 
1067  eta_star = eta_k_star*(1.d0 + (4.d0/5.d0)*eps(m)*g0(m) &
1068  * (1.d0+c_e) ) + (3.d0/5.d0)*gamma_star
1069 
1070  eta0 = 5.0d0*m_pm*dsqrt(th(m)/pi) / (16.d0*d_pm*d_pm)
1071 
1072 ! added Ro_g = 0 for granular flows (no gas).
1073  IF(switch == zero .OR. ro_g_avg == zero) THEN
1074  mu_star = eta0
1075  ELSEIF(th(m) .LT. small_number)THEN
1076  mu_star = zero
1077  ELSE
1078  mu_star = ros_avg(m)*eps(m)*g0(m)*th(m)*eta0 / &
1079  ( ros_avg(m)*eps(m)*g0(m)*th(m) + &
1080  (2.d0*dga*eta0/ros_avg(m)) ) ! Note dgA is ~F_gs/ep_s
1081  ENDIF
1082 
1083 ! shear viscosity in Mth solids phase (add to frictional part)
1084  mu_s = mu_star * eta_star
1085 
1086 ! defining granular pressure (for Jenkins BC)
1087  ps = ros_avg(m)*eps(m)*th(m)*(one+2.d0*(one+c_e)*g0eps_avg)
1088  ! ~ROs_avg(m)*EPS(M)*TH(M)*press_star
1089 
1090 
1091  CASE (gtsh_2012)
1092  chi = g0(m)
1093 
1094  rdissp = one
1095  IF(eps(m) > small_number) rdissp = &
1096  one + 3d0*dsqrt(eps(m)/2d0) + 135d0/64d0*eps(m)*dlog(eps(m)) + &
1097  11.26d0*eps(m)*(one-5.1d0*eps(m)+16.57d0*eps(m)**2-21.77d0* &
1098  eps(m)**3) - eps(m)*chi*dlog(epm)
1099 
1100  re_t = ro_g_avg*d_pm*dsqrt(th(m)) / mu_g_avg
1101  re_g = epg*ro_g_avg*d_pm*vrel/mu_g_avg
1102  rdiss = rdissp + re_t * k_phi(eps(m))
1103  gamaavg = 3d0*pi*mu_g_avg*d_pm*rdiss
1104  mu2_0 = dsqrt(2d0*pi) * chi * (one-c_e**2)
1105  mu4_0 = (4.5d0+c_e**2) * mu2_0
1106  mu4_1 = (6.46875d0+0.9375d0*c_e**2)*mu2_0 + 2d0*dsqrt(2d0*pi)* &
1107  chi*(one+c_e)
1108  a2_gtshw = zero ! for eps = zero
1109 
1110  if((eps(m)> small_number) .and. (th(m) > small_number)) then
1111  zeta_star = 4.5d0*dsqrt(2d0*pi)*(ro_g_avg/ros_avg(m))**2*re_g**2 * &
1112  s_star(eps(m),chi) / (eps(m)*(one-eps(m))**2 * re_t**4)
1113  a2_gtshw = (5d0*mu2_0 - mu4_0) / (mu4_1 - 5d0* &
1114  (19d0/16d0*mu2_0 - 1.5d0*zeta_star))
1115  endif
1116 
1117  eta0 = 0.3125d0/(dsqrt(pi)*d_pm**2)*m_pm*dsqrt(th(m))
1118  nu0 = (96.d0/5.d0)*(eps(m)/d_pm)*dsqrt(th(m)/pi)
1119  nun = 0.25d0*nu0*chi*(3d0-c_e)*(one+c_e) * &
1120  (one+0.7375d0*a2_gtshw)
1121  nuk = nu0*(one+c_e)/3d0*chi*( one+2.0625d0*(one-c_e)+ &
1122  ((947d0-579*c_e)/256d0*a2_gtshw) )
1123  edt_s = 4d0/3d0*dsqrt(pi)*(one-c_e**2)*chi* &
1124  (one+0.1875d0*a2_gtshw)*nu_pm*d_pm**2*dsqrt(th(m))
1125 
1126  if((th(m) > small_number) .and. (eps(m) > small_number)) then
1127  zeta_avg = one/6d0*d_pm*vrel**2*(3d0*pi*mu_g_avg*d_pm/&
1128  m_pm)**2 / dsqrt(th(m)) * s_star(eps(m), chi)
1129  etak = ros_avg(m)*eps(m)*th(m) / (nun-0.5d0*( &
1130  edt_s-zeta_avg/th(m) - 2d0*gamaavg/m_pm)) * &
1131  (one -0.4d0 * (one+c_e)*(one-3d0*c_e)*eps(m)*chi)
1132  else
1133  etak = zero
1134  endif
1135 
1136  mu_b_avg = 25.6d0/pi * eps(m)**2 * chi *(one+c_e) * &
1137  (one - a2_gtshw/16d0)*eta0
1138 
1139  mu_s = etak*(one+0.8d0*eps(m)*chi*(one+c_e)) + &
1140  0.6d0*mu_b_avg
1141 
1142 ! defining granular pressure (for Jenkins BC)
1143  ps = ros_avg(m)*eps(m)*th(m)*(one+2.d0*(one+c_e)*g0eps_avg)
1144 
1145  CASE DEFAULT
1146 ! should never hit this
1147  WRITE (*, '(A)') 'CALC_GRBDRY => F_HW '
1148  WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', kt_type
1149  call mfix_exit(mype)
1150  END SELECT
1151 
1152 ! setting the coefficients for JJ BC
1153 ! -------------------------------------------------------------------
1154  SELECT CASE (kt_type_enum)
1155 
1156  CASE (lun_1984, simonin_1996, ahmadi_1995, gd_1999, gtsh_2012)
1157 ! KT without particle mass in their definition of granular temperature
1158 ! and theta = granular temperature OR
1159 ! KT with particle mass in their definition of granular temperature
1160 ! and theta = granular temperature/mass
1161 
1162 ! Jenkins BC is implemented for these KT
1163  IF(jenkins) THEN
1164  IF (vslip == zero) THEN
1165 ! if solids velocity field is initialized to zero, use free slip bc
1166  f_2 = zero
1167  ELSE
1168 ! As I understand from soil mechanic papers, the coefficient mu in
1169 ! Jenkins paper is tan_Phi_w. T. See for example, G.I. Tardos, PT,
1170 ! 92 (1997), 61-74, equation (1). sof
1171 
1172 ! here F_2 divided by VSLIP to use the same bc as Johnson&Jackson
1173  f_2 = tan_phi_w*ps/vslip
1174 
1175  ENDIF
1176  ELSE
1177  IF(.NOT. bc_jj_m) THEN
1178  f_2 = (phip*dsqrt(3d0*th(m))*pi*ros_avg(m)*eps(m)*&
1179  g0(m))/(6d0*(one-ep_star_avg))
1180  ELSE
1181  f_2 = (phip_jj(vslip,th(m))*dsqrt(3d0*th(m))*pi*&
1182  ros_avg(m)*eps(m)*g0(m))/(6d0*(one-ep_star_avg))
1183  ENDIF
1184  ENDIF ! end if(Jenkins)/else
1185 
1186 
1187  CASE (ia_2005)
1188 ! KTs with particle mass in their definition of granular temperature
1189 ! and theta = granular temperature
1190  IF(.NOT. bc_jj_m) THEN
1191  f_2 = (phip*dsqrt(3.d0*th(m)/m_pm)*pi*ros_avg(m)*eps(m)*&
1192  g0(m))/(6.d0*(one-ep_star_avg))
1193  ELSE
1194  f_2 = (phip_jj(vslip,th(m))*dsqrt(3.d0*th(m)/m_pm)*pi*&
1195  ros_avg(m)*eps(m)*g0(m))/(6.d0*(one-ep_star_avg))
1196  ENDIF
1197 
1198 
1199  CASE DEFAULT
1200 ! should never hit this
1201  WRITE (*, '(A)') 'CALC_GRBDRY => F_HW '
1202  WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', kt_type
1203  call mfix_exit(mype)
1204  END SELECT
1205 
1206  f_hw = f_2/mu_s
1207 
1208  RETURN
1209  END FUNCTION f_hw
1210 
1211 
1212 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1213 ! C
1214 ! Subroutine CG_CALC_GRBDRY C
1215 ! Purpose: Cut cell version C
1216 ! C
1217 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1218 
1219  SUBROUTINE cg_calc_grbdry(IJK, TYPE_OF_CELL, M, L, F_2)
1220 
1221 !-----------------------------------------------
1222 ! Modules
1223 !-----------------------------------------------
1224  USE bc
1225  USE compar
1226  USE constant
1227  USE fldvar
1228  USE fun_avg
1229  USE functions
1230  USE geometry
1231  USE indices
1232  USE param
1233  USE param1
1234  USE physprop
1235  USE rdf
1236  USE run
1237  USE turb
1238  USE visc_s
1239  Use cutcell
1240  use toleranc
1241  IMPLICIT NONE
1242 !-----------------------------------------------
1243 ! Dummy Arguments
1244 !-----------------------------------------------
1245 ! IJK indices
1246  INTEGER, INTENT(IN) :: IJK
1247 
1248  CHARACTER (LEN=*) :: TYPE_OF_CELL
1249 ! Solids phase index
1250  INTEGER, INTENT(IN) :: M
1251 ! Index corresponding to boundary condition
1252  INTEGER, INTENT(IN) :: L
1253 
1254  DOUBLE PRECISION :: F_2
1255 
1256 !-----------------------------------------------
1257 ! Local Variables
1258 !-----------------------------------------------
1259 ! IJK indices
1260  INTEGER IJK1, IJK2
1261 ! Other indices
1262  INTEGER IJK2E, IPJK2, IPJMK2
1263 !
1264  DOUBLE PRECISION :: smallTheta
1265 ! Average scalars
1266  DOUBLE PRECISION EPg_avg, Mu_g_avg, RO_g_avg
1267 ! Average void fraction at packing
1268  DOUBLE PRECISION ep_star_avg
1269 ! Average scalars modified to include all solid phases
1270  DOUBLE PRECISION EPs_avg(dimension_m), DP_avg(dimension_m),&
1271  TH_avg(dimension_m)
1272 ! Average solids density
1273  DOUBLE PRECISION ROS_AVG(dimension_m)
1274 ! Average Simonin and Ahmadi variables (sof)
1275  DOUBLE PRECISION K_12_avg, Tau_12_avg, Tau_1_avg
1276 ! Average velocities
1277  DOUBLE PRECISION WGC1, WGC2
1278 ! Solids phase index
1279  INTEGER MM
1280 ! values of U_sm, V_sm, W_sm at appropriate place on boundary wall
1281  DOUBLE PRECISION USCM, VSCM,WSCM
1282  DOUBLE PRECISION WSCM1,WSCM2
1283 ! values of U_g, V_g, W_g at appropriate place on boundary wall
1284  DOUBLE PRECISION UGC, VGC, WGC
1285 ! Magnitude of gas-solids relative velocity
1286  DOUBLE PRECISION VREL
1287 ! slip velocity between wall and particles for Jenkins bc (sof)
1288  DOUBLE PRECISION VSLIP
1289 ! radial distribution function at contact
1290  DOUBLE PRECISION g0(dimension_m)
1291 ! Sum of eps*G_0
1292  DOUBLE PRECISION g0EPs_avg
1293 !-----------------------------------------------
1294 
1295 ! Note: EP_s, MU_g, and RO_g are undefined at IJK1 (wall cell). Hence
1296 ! IJK2 (fluid cell) is used in averages.
1297 !
1298 
1299  smalltheta = (to_si)**4 * zero_ep_s
1300 
1301  SELECT CASE (type_of_cell)
1302 
1303  CASE('U_MOMENTUM')
1304 
1305  ijk2 = east_of(ijk)
1306  ijk2e = east_of(ijk2)
1307 
1308  epg_avg = (vol(ijk)*ep_g(ijk) + vol(ijk2)*ep_g(ijk2))/(vol(ijk) + vol(ijk2))
1309  ep_star_avg = (vol(ijk)*ep_star_array(ijk) + vol(ijk2)*ep_star_array(ijk2))/(vol(ijk) + vol(ijk2))
1310  mu_g_avg = (vol(ijk)*mu_g(ijk) + vol(ijk2)*mu_g(ijk2))/(vol(ijk) + vol(ijk2))
1311  ro_g_avg = (vol(ijk)*ro_g(ijk) + vol(ijk2)*ro_g(ijk2))/(vol(ijk) + vol(ijk2))
1312 
1313  k_12_avg = zero
1314  tau_12_avg = zero
1315  tau_1_avg = zero
1316  IF(kt_type_enum == simonin_1996 .OR.&
1317  kt_type_enum == ahmadi_1995) THEN ! not converted to CG
1318  k_12_avg = avg_x(k_12(ijk2), k_12(ijk2e), i_of(ijk2))
1319  tau_12_avg = avg_x(tau_12(ijk2), tau_12(ijk2e), i_of(ijk2))
1320  tau_1_avg = avg_x(tau_1(ijk2), tau_1(ijk2e), i_of(ijk2))
1321  ENDIF
1322 
1323  g0eps_avg = zero
1324  DO mm = 1, smax
1325  g0(mm) = g_0avg(ijk, ijk, 'X', i_of(ijk), m, mm)
1326  eps_avg(mm) = (vol(ijk)*ep_s(ijk, mm) + vol(ijk2)*ep_s(ijk2, mm))/(vol(ijk) + vol(ijk2))
1327  dp_avg(mm) = (vol(ijk)*d_p(ijk, mm) + vol(ijk2)*d_p(ijk2, mm))/(vol(ijk) + vol(ijk2))
1328  g0eps_avg = g0eps_avg + g_0avg(ijk, ijk, 'X', i_of(ijk), m, mm) &
1329  * (vol(ijk)*ep_s(ijk, mm) + vol(ijk2)*ep_s(ijk2, mm))/(vol(ijk) + vol(ijk2))
1330  ros_avg(mm) = (vol(ijk)*ro_s(ijk, mm) + vol(ijk2)*ro_s(ijk2, mm))/(vol(ijk) + vol(ijk2))
1331 
1332  th_avg(mm) = (vol(ijk)*theta_m(ijk,mm) + vol(ijk2)*theta_m(ijk2,mm))/(vol(ijk) + vol(ijk2))
1333  IF(th_avg(mm) < zero) th_avg(mm) = smalltheta
1334  ENDDO
1335 
1336 
1337 ! Calculate velocity components at i+1/2, j+1/2, k (relative to IJK1) ! not converted to CG
1338  ugc = avg_y(u_g(ijk1), u_g(ijk2),j_of(ijk1))
1339  vgc = avg_x(v_g(ijk1), v_g(ipjmk2),i_of(ijk1))
1340  wgc1 = avg_x(avg_z_t(w_g(km_of(ijk2)), w_g(ijk2)),&
1341  avg_z_t(w_g(km_of(ipjk2)), w_g(ipjk2)),&
1342  i_of(ijk2))
1343  wgc2 = avg_x(avg_z_t(w_g(km_of(ijk1)), w_g(ijk1)),&
1344  avg_z_t(w_g(km_of(ipjmk2)), w_g(ipjmk2)),&
1345  i_of(ijk1))
1346  wgc = avg_y(wgc2, wgc1, j_of(ijk1))
1347  uscm = avg_y(u_s(ijk1,m), u_s(ijk2,m),j_of(ijk1))
1348  vscm = avg_x(v_s(ijk1,m), v_s(ipjmk2,m),i_of(ijk1))
1349  wscm1= avg_x(avg_z_t(w_s(km_of(ijk2),m),w_s(ijk2,m)),&
1350  avg_z_t(w_s(km_of(ipjk2),m),w_s(ipjk2,m)),&
1351  i_of(ijk2))
1352  wscm2= avg_x(avg_z_t(w_s(km_of(ijk1),m),w_s(ijk1,m)),&
1353  avg_z_t(w_s(km_of(ipjmk2),m),w_s(ipjmk2,m)),&
1354  i_of(ijk1))
1355  wscm = avg_y(wscm2, wscm1, j_of(ijk1))
1356 
1357 ! magnitude of gas-solids relative velocity
1358  vrel = dsqrt( (ugc - uscm)**2 + (vgc - vscm)**2 + &
1359  (wgc - wscm)**2 )
1360 
1361 ! slip velocity for use in Jenkins bc (sof)
1362  vslip= dsqrt( (uscm-bc_uw_s(l,m))**2 + (vscm-bc_vw_s(l,m))**2 &
1363  + (wscm-bc_ww_s(l,m))**2 )
1364 
1365  CALL get_cg_f2(g0, eps_avg, epg_avg, ep_star_avg, &
1366  g0eps_avg, th_avg, mu_g_avg, ro_g_avg, ros_avg,&
1367  dp_avg, k_12_avg, tau_12_avg, tau_1_avg, &
1368  vrel, vslip, m,f_2)
1369 
1370 
1371  CASE('V_MOMENTUM')
1372 
1373  ijk2 = north_of(ijk)
1374 
1375  epg_avg = (vol(ijk)*ep_g(ijk) + vol(ijk2)*ep_g(ijk2))/(vol(ijk) + vol(ijk2))
1376  ep_star_avg = (vol(ijk)*ep_star_array(ijk) + vol(ijk2)*ep_star_array(ijk2))/(vol(ijk) + vol(ijk2))
1377  mu_g_avg = (vol(ijk)*mu_g(ijk) + vol(ijk2)*mu_g(ijk2))/(vol(ijk) + vol(ijk2))
1378  ro_g_avg = (vol(ijk)*ro_g(ijk) + vol(ijk2)*ro_g(ijk2))/(vol(ijk) + vol(ijk2))
1379 
1380  k_12_avg = zero
1381  tau_12_avg = zero
1382  tau_1_avg = zero
1383  IF(kt_type_enum == simonin_1996 .OR.&
1384  kt_type_enum == ahmadi_1995) THEN ! not converted to CG
1385  tau_1_avg = avg_x(tau_1(ijk2), tau_1(ijk2e), i_of(ijk2))
1386 ! Simonin only:
1387  k_12_avg = avg_x(k_12(ijk2), k_12(ijk2e), i_of(ijk2))
1388  tau_12_avg = avg_x(tau_12(ijk2), tau_12(ijk2e), i_of(ijk2))
1389  ENDIF
1390 
1391  g0eps_avg = zero
1392  DO mm = 1, smax
1393  g0(mm) = g_0avg(ijk, ijk, 'X', i_of(ijk), m, mm)
1394  eps_avg(mm) = (vol(ijk)*ep_s(ijk, mm) + vol(ijk2)*ep_s(ijk2, mm))/(vol(ijk) + vol(ijk2))
1395  dp_avg(mm) = (vol(ijk)*d_p(ijk, mm) + vol(ijk2)*d_p(ijk2, mm))/(vol(ijk) + vol(ijk2))
1396  g0eps_avg = g0eps_avg + g_0avg(ijk, ijk, 'X', i_of(ijk), m, mm) &
1397  * (vol(ijk)*ep_s(ijk, mm) + vol(ijk2)*ep_s(ijk2, mm))/(vol(ijk) + vol(ijk2))
1398  ros_avg(mm) = (vol(ijk)*ro_s(ijk, mm) + vol(ijk2)*ro_s(ijk2, mm))/(vol(ijk) + vol(ijk2))
1399  th_avg(mm) = (vol(ijk)*theta_m(ijk,mm) + vol(ijk2)*theta_m(ijk2,mm))/(vol(ijk) + vol(ijk2))
1400  IF(th_avg(mm) < zero) th_avg(mm) = smalltheta
1401  ENDDO
1402 
1403 ! Calculate velocity components at i+1/2, j+1/2, k (relative to IJK1) ! not converted to CG
1404  ugc = avg_y(u_g(ijk1), u_g(ijk2),j_of(ijk1))
1405  vgc = avg_x(v_g(ijk1), v_g(ipjmk2),i_of(ijk1))
1406  wgc1 = avg_x(avg_z_t(w_g(km_of(ijk2)), w_g(ijk2)),&
1407  avg_z_t(w_g(km_of(ipjk2)), w_g(ipjk2)),&
1408  i_of(ijk2))
1409  wgc2 = avg_x(avg_z_t(w_g(km_of(ijk1)), w_g(ijk1)),&
1410  avg_z_t(w_g(km_of(ipjmk2)), w_g(ipjmk2)),&
1411  i_of(ijk1))
1412  wgc = avg_y(wgc2, wgc1, j_of(ijk1))
1413  uscm = avg_y(u_s(ijk1,m), u_s(ijk2,m),j_of(ijk1))
1414  vscm = avg_x(v_s(ijk1,m), v_s(ipjmk2,m),i_of(ijk1))
1415  wscm1= avg_x(avg_z_t(w_s(km_of(ijk2),m),w_s(ijk2,m)),&
1416  avg_z_t(w_s(km_of(ipjk2),m),w_s(ipjk2,m)),&
1417  i_of(ijk2))
1418  wscm2= avg_x(avg_z_t(w_s(km_of(ijk1),m),w_s(ijk1,m)),&
1419  avg_z_t(w_s(km_of(ipjmk2),m),w_s(ipjmk2,m)),&
1420  i_of(ijk1))
1421  wscm = avg_y(wscm2, wscm1, j_of(ijk1))
1422 
1423 ! magnitude of gas-solids relative velocity
1424 
1425  vrel = dsqrt( (ugc - uscm)**2 + (vgc - vscm)**2 + &
1426  (wgc - wscm)**2 )
1427 
1428 ! slip velocity for use in Jenkins bc (sof)
1429  vslip= dsqrt( (uscm-bc_uw_s(l,m))**2 + (vscm-bc_vw_s(l,m))**2 &
1430  + (wscm-bc_ww_s(l,m))**2 )
1431 
1432 
1433  CALL get_cg_f2(g0, eps_avg, epg_avg, ep_star_avg, &
1434  g0eps_avg, th_avg, mu_g_avg, ro_g_avg, ros_avg,&
1435  dp_avg, k_12_avg, tau_12_avg, tau_1_avg, &
1436  vrel, vslip, m,f_2)
1437 
1438 
1439  CASE('W_MOMENTUM')
1440 
1441  ijk2 = top_of(ijk)
1442 
1443  epg_avg = (vol(ijk)*ep_g(ijk) + vol(ijk2)*ep_g(ijk2))/(vol(ijk) + vol(ijk2))
1444  ep_star_avg = (vol(ijk)*ep_star_array(ijk) + vol(ijk2)*ep_star_array(ijk2))/(vol(ijk) + vol(ijk2))
1445  mu_g_avg = (vol(ijk)*mu_g(ijk) + vol(ijk2)*mu_g(ijk2))/(vol(ijk) + vol(ijk2))
1446  ro_g_avg = (vol(ijk)*ro_g(ijk) + vol(ijk2)*ro_g(ijk2))/(vol(ijk) + vol(ijk2))
1447 
1448  k_12_avg = zero
1449  tau_12_avg = zero
1450  tau_1_avg = zero
1451  IF(kt_type_enum == simonin_1996 .OR.&
1452  kt_type_enum == ahmadi_1995) THEN ! not converted to CG
1453  tau_1_avg = avg_x(tau_1(ijk2), tau_1(ijk2e), i_of(ijk2))
1454 ! Simonon only:
1455  k_12_avg = avg_x(k_12(ijk2), k_12(ijk2e), i_of(ijk2))
1456  tau_12_avg = avg_x(tau_12(ijk2), tau_12(ijk2e), i_of(ijk2))
1457  ENDIF
1458 
1459  g0eps_avg = zero
1460  DO mm = 1, smax
1461  g0(mm) = g_0avg(ijk, ijk, 'X', i_of(ijk), m, mm)
1462  eps_avg(mm) = (vol(ijk)*ep_s(ijk, mm) + vol(ijk2)*ep_s(ijk2, mm))/(vol(ijk) + vol(ijk2))
1463  dp_avg(mm) = (vol(ijk)*d_p(ijk, mm) + vol(ijk2)*d_p(ijk2, mm))/(vol(ijk) + vol(ijk2))
1464  ros_avg(mm) = (vol(ijk)*ro_s(ijk, mm) + vol(ijk2)*ro_s(ijk2, mm))/(vol(ijk) + vol(ijk2))
1465  g0eps_avg = g0eps_avg + g_0avg(ijk, ijk, 'X', i_of(ijk), m, mm) &
1466  * (vol(ijk)*ep_s(ijk, mm) + vol(ijk2)*ep_s(ijk2, mm))/(vol(ijk) + vol(ijk2))
1467 
1468  th_avg(mm) = (vol(ijk)*theta_m(ijk,mm) + vol(ijk2)*theta_m(ijk2,mm))/(vol(ijk) + vol(ijk2))
1469  IF(th_avg(mm) < zero) th_avg(mm) = smalltheta
1470  ENDDO
1471 
1472 ! Calculate velocity components at i+1/2, j+1/2, k (relative to IJK1) ! not converted to CG
1473  ugc = avg_y(u_g(ijk1), u_g(ijk2),j_of(ijk1))
1474  vgc = avg_x(v_g(ijk1), v_g(ipjmk2),i_of(ijk1))
1475  wgc1 = avg_x(avg_z_t(w_g(km_of(ijk2)), w_g(ijk2)),&
1476  avg_z_t(w_g(km_of(ipjk2)), w_g(ipjk2)),&
1477  i_of(ijk2))
1478  wgc2 = avg_x(avg_z_t(w_g(km_of(ijk1)), w_g(ijk1)),&
1479  avg_z_t(w_g(km_of(ipjmk2)), w_g(ipjmk2)),&
1480  i_of(ijk1))
1481  wgc = avg_y(wgc2, wgc1, j_of(ijk1))
1482  uscm = avg_y(u_s(ijk1,m), u_s(ijk2,m),j_of(ijk1))
1483  vscm = avg_x(v_s(ijk1,m), v_s(ipjmk2,m),i_of(ijk1))
1484  wscm1= avg_x(avg_z_t(w_s(km_of(ijk2),m),w_s(ijk2,m)),&
1485  avg_z_t(w_s(km_of(ipjk2),m),w_s(ipjk2,m)),&
1486  i_of(ijk2))
1487  wscm2= avg_x(avg_z_t(w_s(km_of(ijk1),m),w_s(ijk1,m)),&
1488  avg_z_t(w_s(km_of(ipjmk2),m),w_s(ipjmk2,m)),&
1489  i_of(ijk1))
1490  wscm = avg_y(wscm2, wscm1, j_of(ijk1))
1491 
1492 ! magnitude of gas-solids relative velocity
1493  vrel = dsqrt( (ugc - uscm)**2 + (vgc - vscm)**2 + &
1494  (wgc - wscm)**2 )
1495 
1496 ! slip velocity for use in Jenkins bc (sof)
1497  vslip= dsqrt( (uscm-bc_uw_s(l,m))**2 + (vscm-bc_vw_s(l,m))**2 &
1498  + (wscm-bc_ww_s(l,m))**2 )
1499 
1500 
1501 ! why bother with this since it should not come here...
1502  CALL get_cg_f2(g0, eps_avg, epg_avg, ep_star_avg, &
1503  g0eps_avg, th_avg, mu_g_avg, ro_g_avg, ros_avg,&
1504  dp_avg, k_12_avg, tau_12_avg, tau_1_avg, &
1505  vrel, vslip, m,f_2)
1506 
1507 
1508  CASE DEFAULT
1509  WRITE(*,*) 'CG_CALC_GRBDRY'
1510  WRITE(*,*) 'UNKNOWN TYPE OF CELL:',type_of_cell
1511  WRITE(*,*) 'ACCEPTABLE TYPES ARE:'
1512  WRITE(*,*) 'U_MOMENTUM'
1513  WRITE(*,*) 'V_MOMENTUM'
1514  WRITE(*,*) 'W_MOMENTUM'
1515  CALL mfix_exit(mype)
1516  END SELECT
1517 
1518 
1519  RETURN
1520 
1521  END SUBROUTINE cg_calc_grbdry
1522 
1523 
1524 
1525 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1526 ! C
1527 ! Subroutine: GET_CG_F2 C
1528 ! Purpose: Compute F_2 for cut cell version C
1529 ! C
1530 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1531  SUBROUTINE get_cg_f2(g0,EPS,EPG, ep_star_avg, &
1532  g0eps_avg,th,mu_g_avg,ro_g_avg,ros_avg, &
1533  dp_avg,k_12_avg, tau_12_avg, tau_1_avg, &
1534  vrel, vslip, m,f_2)
1535 
1536 !-----------------------------------------------
1537 ! Modules
1538 !-----------------------------------------------
1539  USE param
1540  USE param1
1541  USE constant
1542  USE physprop
1543  USE run
1544  USE fldvar
1545  USE mpi_utility
1546  Use cutcell
1547  IMPLICIT NONE
1548 !-----------------------------------------------
1549 ! Dummy Arguments
1550 !-----------------------------------------------
1551 ! Radial distribution function of solids phase M with each
1552 ! other solids phase
1553  DOUBLE PRECISION, INTENT(IN) :: g0(dimension_m)
1554 ! Average solids volume fraction of each solids phase
1555  DOUBLE PRECISION, INTENT(IN) :: EPS(dimension_m)
1556 ! Average solids and gas volume fraction
1557  DOUBLE PRECISION, INTENT(IN) :: EPG, ep_star_avg
1558 ! Sum of eps*G_0
1559  DOUBLE PRECISION, INTENT(INOUT) :: g0EPs_avg
1560 ! Average theta_m
1561  DOUBLE PRECISION, INTENT(INOUT) :: TH (dimension_m)
1562 ! Average gas viscosity
1563  DOUBLE PRECISION, INTENT(IN) :: Mu_g_avg
1564 ! Average gas density
1565  DOUBLE PRECISION, INTENT(IN) :: RO_g_avg
1566 !QX: Average solids density
1567  DOUBLE PRECISION, INTENT(IN) :: ROS_avg(dimension_m)
1568 ! Average particle diameter of each solids phase
1569  DOUBLE PRECISION, INTENT(IN) :: DP_avg(dimension_m)
1570 ! Average cross-correlation K_12 and lagrangian integral time-scale
1571  DOUBLE PRECISION, INTENT(IN) :: K_12_avg, Tau_12_avg, Tau_1_avg
1572 ! Magnitude of slip velocity between two phases
1573  DOUBLE PRECISION, INTENT(IN) :: VREL
1574 ! Slip velocity between wall and particles
1575  DOUBLE PRECISION, INTENT(IN) :: VSLIP
1576 ! Solids phase index
1577  INTEGER, INTENT(IN) :: M
1578 !-----------------------------------------------
1579 ! Local Variables
1580 !-----------------------------------------------
1581 ! Solids pressure
1582  DOUBLE PRECISION :: Ps
1583  DOUBLE PRECISION :: F_2
1584  DOUBLE PRECISION :: D_PM, M_PM
1585 !-----------------------------------------------
1586 ! Functions
1587 !-----------------------------------------------
1588 ! Variable specularity coefficient
1589  DOUBLE PRECISION :: PHIP_JJ
1590 !-----------------------------------------------
1591 
1592 ! This is done here similar to bc_theta to avoid small negative values of
1593 ! Theta coming most probably from linear solver
1594  IF(th(m) .LE. zero)THEN
1595  th(m) = 1d-8
1596  IF (mype.eq.pe_io) THEN
1597  WRITE(*,*) &
1598  'Warning: Negative granular temp at wall set to 1e-8'
1599 ! CALL WRITE_ERROR('THETA_HW_CW', LINE, 1)
1600  ENDIF
1601  ENDIF
1602 
1603 ! common variables
1604  d_pm = dp_avg(m)
1605  m_pm = (pi/6.d0)*(d_pm**3)*ros_avg(m)
1606 
1607 ! Defining viscosity according to each KT for later use in JJ BC
1608 ! Also defining solids pressure according to each KT for use if JENKINS
1609 ! -------------------------------------------------------------------
1610  SELECT CASE (kt_type_enum)
1611  CASE (lun_1984)
1612 ! defining granular pressure (for Jenkins BC)
1613  ps = ros_avg(m)*eps(m)*th(m)*(one+4.d0*eta*g0eps_avg)
1614 
1615 
1616  CASE (simonin_1996)
1617 ! defining granular pressure (for Jenkins BC)
1618  ps = ros_avg(m)*eps(m)*th(m)*(one+4.d0*eta*g0eps_avg)
1619 
1620 
1621  CASE (ahmadi_1995)
1622 ! defining granular pressure (for Jenkins BC)
1623  ps = ros_avg(m)*eps(m)*th(m)*&
1624  ((one + 4.0d0*g0eps_avg)+half*(one -c_e*c_e))
1625 
1626 
1627  CASE(ia_2005)
1628 ! Use original IA theory if SWITCH_IA is false
1629 ! no ps since no jenkins for this theory
1630 
1631  CASE(gd_1999)
1632 
1633 ! defining granular pressure (for Jenkins BC)
1634  ps = ros_avg(m)*eps(m)*th(m)*(one+2.d0*(one+c_e)*g0eps_avg)
1635  ! ~ROs_avg(m)*EPS(M)*TH(M)*press_star
1636 
1637 
1638  CASE (gtsh_2012)
1639 ! defining granular pressure (for Jenkins BC)
1640  ps = ros_avg(m)*eps(m)*th(m)*(one+2.d0*(one+c_e)*g0eps_avg)
1641 
1642 
1643  CASE DEFAULT
1644 ! should never hit this
1645  WRITE (*, '(A)') 'CG_CALC_GRBDRY => GET_CG_F2 '
1646  WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', kt_type
1647  call mfix_exit(mype)
1648  END SELECT
1649 
1650 
1651 ! setting the coefficients for JJ BC
1652 ! -------------------------------------------------------------------
1653  SELECT CASE (kt_type_enum)
1654 
1655  CASE (lun_1984, simonin_1996, ahmadi_1995, gd_1999, gtsh_2012)
1656 ! KT without particle mass in their definition of granular temperature
1657 ! and theta = granular temperature OR
1658 ! KT with particle mass in their definition of granular temperature
1659 ! and theta = granular temperature/mass
1660 
1661 ! Jenkins BC is implemented for these KT
1662  IF(jenkins) THEN
1663  IF (vslip == zero) THEN
1664 ! if solids velocity field is initialized to zero, use free slip bc
1665  f_2 = zero
1666  ELSE
1667 ! As I understand from soil mechanic papers, the coefficient mu in
1668 ! Jenkins paper is tan_Phi_w. T. See for example, G.I. Tardos, PT,
1669 ! 92 (1997), 61-74, equation (1). sof
1670 
1671 ! here F_2 divided by VSLIP to use the same bc as Johnson&Jackson
1672  f_2 = tan_phi_w*ps/vslip
1673 
1674  ENDIF
1675  ELSE
1676  IF(.NOT. bc_jj_m) THEN
1677  f_2 = (phip*dsqrt(3d0*th(m))*pi*ros_avg(m)*eps(m)*&
1678  g0(m))/(6d0*(one-ep_star_avg))
1679  ELSE
1680  f_2 = (phip_jj(vslip,th(m))*dsqrt(3d0*th(m))*pi*&
1681  ros_avg(m)*eps(m)*g0(m))/(6d0*(one-ep_star_avg))
1682  ENDIF
1683  ENDIF ! end if(Jenkins)/else
1684 
1685  CASE(ia_2005)
1686 ! KTs with particle mass in their definition of granular temperature
1687 ! and theta = granular temperature
1688  IF(.NOT. bc_jj_m) THEN
1689  f_2 = (phip*dsqrt(3.d0*th(m)/m_pm)*pi*ros_avg(m)*&
1690  eps(m)*g0(m))/(6.d0*(one-ep_star_avg))
1691  ELSE
1692  f_2 = (phip_jj(vslip,th(m))*dsqrt(3.d0*th(m)/m_pm)*&
1693  pi*ros_avg(m)*eps(m)*g0(m))/(6.d0*(one-ep_star_avg))
1694  ENDIF
1695 
1696 
1697  CASE DEFAULT
1698 ! should never hit this
1699  WRITE (*, '(A)') 'CALC_GRBDRY => F_HW '
1700  WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', kt_type
1701  call mfix_exit(mype)
1702  END SELECT
1703 
1704 ! all the code used to calculate Mu_s was deleted from this routine
1705 ! F_HW = F_2/Mu_s ! Only F_2 is actually needed
1706 
1707 
1708  RETURN
1709  END SUBROUTINE get_cg_f2
1710 END MODULE calc_gr_boundary
bc::bc_ww_s
double precision, dimension(dimension_bc, dim_m) bc_ww_s
Definition: bc_mod.f:328
mpi_utility
Definition: mpi_utility_mod.f:6
constant::c_e
double precision c_e
Definition: constant_mod.f:105
fldvar::v_s
double precision, dimension(:,:), allocatable v_s
Definition: fldvar_mod.f:105
constant::to_si
double precision to_si
Definition: constant_mod.f:146
turb::tau_1
double precision, dimension(:), allocatable tau_1
Definition: turb_mod.f:19
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
visc_s
Definition: visc_s_mod.f:1
bc::bc_uw_s
double precision, dimension(dimension_bc, dim_m) bc_uw_s
Definition: bc_mod.f:322
fldvar::ep_g
double precision, dimension(:), allocatable ep_g
Definition: fldvar_mod.f:15
param1::one
double precision, parameter one
Definition: param1_mod.f:29
write_error
subroutine write_error(NAME, LINE, LMAX)
Definition: write_error.f:21
fldvar::w_s
double precision, dimension(:,:), allocatable w_s
Definition: fldvar_mod.f:117
run::friction
logical friction
Definition: run_mod.f:149
constant::switch
double precision, parameter switch
Definition: constant_mod.f:53
calc_gr_boundary::f_hw
double precision function f_hw(g0, EPS, EPG, ep_star_avg, g0EPs_avg, TH, Mu_g_avg, RO_g_avg, ROs_avg, DP_avg, K_12_avg, Tau_12_avg, Tau_1_avg, VREL, VSLIP, M)
Definition: calc_grbdry.f:794
turb::k_12
double precision, dimension(:), allocatable k_12
Definition: turb_mod.f:17
kintheory::s_star
double precision function s_star(phi, Chi)
Definition: kintheory_mod.f:822
turb::tau_12
double precision, dimension(:), allocatable tau_12
Definition: turb_mod.f:18
fldvar::u_s
double precision, dimension(:,:), allocatable u_s
Definition: fldvar_mod.f:93
calc_gw_hw_cw
subroutine calc_gw_hw_cw(g0, EPS, EPG, ep_star_avg, g0EP_avg, TH, Mu_g_avg, RO_g_avg, ROs_avg, DP_avg, K_12_avg, Tau_12_avg, Tau_1_avg, VREL, VSLIP, DEL_U, S_S, S_dd, VEL, W_VEL, M, gw, hw, cw)
Definition: calc_u_friction.f:32
visc_s::ep_star_array
double precision, dimension(:), allocatable ep_star_array
Definition: visc_s_mod.f:54
turb
Definition: turb_mod.f:2
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
fun_avg
Definition: fun_avg_mod.f:1
calc_gr_boundary::calc_grbdry
subroutine calc_grbdry(IJK1, IJK2, FCELL, COM, M, L, Gw, Hw, Cw)
Definition: calc_grbdry.f:25
run::jenkins
logical jenkins
Definition: run_mod.f:166
fldvar::d_p
double precision, dimension(:,:), allocatable d_p
Definition: fldvar_mod.f:57
constant::alpha
double precision, parameter alpha
Definition: constant_mod.f:62
calc_s_ddot_s
subroutine calc_s_ddot_s(IJK1, IJK2, FCELL, COM, M, DEL_DOT_U, S_DDOT_S,
Definition: calc_s_ddot_s.f:24
exit::mfix_exit
subroutine mfix_exit(myID, normal_termination)
Definition: exit.f:5
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
run::bc_jj_m
logical bc_jj_m
Definition: run_mod.f:168
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
exit
Definition: exit.f:2
constant::eta
double precision eta
Definition: constant_mod.f:108
toleranc::zero_ep_s
double precision, parameter zero_ep_s
Definition: toleranc_mod.f:15
fldvar::theta_m
double precision, dimension(:,:), allocatable theta_m
Definition: fldvar_mod.f:149
constant::eps_f_min
double precision eps_f_min
Definition: constant_mod.f:100
fldvar::v_g
double precision, dimension(:), allocatable v_g
Definition: fldvar_mod.f:99
constant::phip
double precision phip
Definition: constant_mod.f:79
fldvar::w_g
double precision, dimension(:), allocatable w_g
Definition: fldvar_mod.f:111
param1::half
double precision, parameter half
Definition: param1_mod.f:28
run
Definition: run_mod.f:13
param1::large_number
double precision, parameter large_number
Definition: param1_mod.f:23
param
Definition: param_mod.f:2
fldvar::ro_s
double precision, dimension(:,:), allocatable ro_s
Definition: fldvar_mod.f:45
kintheory::k_phi
double precision function k_phi(phi)
Definition: kintheory_mod.f:777
compar::mype
integer mype
Definition: compar_mod.f:24
physprop::mu_g
double precision, dimension(:), allocatable mu_g
Definition: physprop_mod.f:68
kintheory::epm
double precision, parameter epm
Definition: kintheory_mod.f:81
fldvar::u_g
double precision, dimension(:), allocatable u_g
Definition: fldvar_mod.f:87
fldvar::ep_s
double precision function ep_s(IJK, xxM)
Definition: fldvar_mod.f:178
constant::tan_phi_w
double precision tan_phi_w
Definition: constant_mod.f:132
ps
Definition: ps_mod.f:22
rdf
Definition: rdf_mod.f:11
physprop::smax
integer smax
Definition: physprop_mod.f:22
calc_gr_boundary::get_cg_f2
subroutine get_cg_f2(g0, EPS, EPG, ep_star_avg, g0EPs_avg, TH, Mu_g_avg, RO_g_avg, Ros_avg, DP_avg, K_12_avg, Tau_12_avg, Tau_1_avg, VREL, VSLIP, M, F_2)
Definition: calc_grbdry.f:1535
calc_gr_boundary
Definition: calc_grbdry.f:1
mpi_utility::exitmpi
subroutine exitmpi(myid)
Definition: mpi_utility_mod.f:5437
constant::pi
double precision, parameter pi
Definition: constant_mod.f:158
geometry::vol
double precision, dimension(:), allocatable vol
Definition: geometry_mod.f:212
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
calc_gr_boundary::cg_calc_grbdry
subroutine cg_calc_grbdry(IJK, TYPE_OF_CELL, M, L, F_2)
Definition: calc_grbdry.f:1220
fldvar::ro_g
double precision, dimension(:), allocatable ro_g
Definition: fldvar_mod.f:32
rdf::g_0avg
double precision function g_0avg(IJK1, IJK2, DIR, L, M1, M2)
Definition: rdf_mod.f:21
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
bc
Definition: bc_mod.f:23
constant::switch_ia
logical, parameter switch_ia
Definition: constant_mod.f:74