MFIX  2016-1
bc_theta.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: BC_THETA C
4 ! C
5 ! Purpose: Calculate boundary conditions for the temperature/ C
6 ! granular energy equation C
7 ! C
8 ! Author: Kapil Agrawal, Princeton University Date: 14-MAR-98 C
9 ! Reviewer: Date: C
10 ! C
11 ! C
12 ! Literature/Document References: C
13 ! Johnson, P. C., and Jackson, R., Frictional-collisional C
14 ! constitutive relations for granular materials, with C
15 ! application to plane shearing, Journal of Fluid Mechanics, C
16 ! 1987, 176, pp. 67-93 C
17 ! C
18 ! Notes: C
19 ! If BC_JJ_PS = 3, then set the wall temperature equal to the C
20 ! adjacent fluid cell temperature (i.e. zero flux) C
21 ! IF BC_JJ_PS = 1, then set the wall temperature based on C
22 ! Johnson and Jackson Pseudo-thermal temp BC C
23 ! IF BC_JJ_PS !=1, !=3 and >0, then set the wall temperature C
24 ! based on Johnson and Jackson Pseudo-thermal temp BC without C
25 ! the generation term C
26 ! C
27 ! C
28 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
29 
30  SUBROUTINE bc_theta(M, A_m, B_m)
31 
32 !-----------------------------------------------
33 ! Modules
34 !-----------------------------------------------
35  USE param
36  USE param1
37  USE parallel
38  USE scales
39  USE constant
40  USE toleranc
41  USE run
42  USE physprop
43  USE fldvar
44  USE visc_s
45  USE geometry
46  USE output
47  USE indices
48  USE bc
49  USE compar
50  USE mpi_utility
51  USE fun_avg
52  USE functions
53  IMPLICIT NONE
54 !-----------------------------------------------
55 ! Dummy arguments
56 !-----------------------------------------------
57 ! Solids phase index
58  INTEGER :: M
59 ! Septadiagonal matrix A_m
60  DOUBLE PRECISION A_m(dimension_3, -3:3, 0:dimension_m)
61 ! Vector b_m
62  DOUBLE PRECISION B_m(dimension_3, 0:dimension_m)
63 !-----------------------------------------------
64 ! Local variables
65 !-----------------------------------------------
66 ! Boundary condition
67  INTEGER :: L
68 ! Indices
69  INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK,&
70  IM, JM, KM
71 ! Granular energy coefficient
72  DOUBLE PRECISION :: Gw, Hw, Cw
73 !-----------------------------------------------
74 
75 ! Setup Johnson and Jackson Pseudo-thermal temp B.C.
76  DO l = 1, dimension_bc
77  IF( bc_defined(l) ) THEN
78  IF(bc_type_enum(l) .EQ. no_slip_wall .OR.&
79  bc_type_enum(l) .EQ. free_slip_wall .OR.&
80  bc_type_enum(l) .EQ. par_slip_wall ) THEN
81  i1 = bc_i_w(l)
82  i2 = bc_i_e(l)
83  j1 = bc_j_s(l)
84  j2 = bc_j_n(l)
85  k1 = bc_k_b(l)
86  k2 = bc_k_t(l)
87 
88  IF(bc_jj_ps(l).GT.0) THEN
89  DO k = k1, k2
90  DO j = j1, j2
91  DO i = i1, i2
92 
93  IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
94  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
95  ijk = funijk(i, j, k)
96  IF (flow_at(ijk)) cycle !checks for pressure outlets
97  im = im1(i)
98  jm = jm1(j)
99  km = km1(k)
100 ! Setting the temperature to zero - recall the center coefficient
101 ! and source vector are negative. The off-diagonal coefficients
102 ! are positive.
103  a_m(ijk, east, m) = zero
104  a_m(ijk, west, m) = zero
105  a_m(ijk, north, m) = zero
106  a_m(ijk, south, m) = zero
107  a_m(ijk, top, m) = zero
108  a_m(ijk, bottom, m) = zero
109  a_m(ijk, 0, m) = -one
110  b_m(ijk, m) = zero
111 
112 ! Checking if the west wall
113  IF(fluid_at(east_of(ijk)))THEN
114  IF(ep_s(east_of(ijk),m).LE.dil_ep_s)THEN
115  a_m(ijk, east, m) = one
116  ELSE
117  IF (bc_jj_ps(l).EQ.3) THEN
118 ! Setting the wall temperature equal to the adjacent fluid cell
119 ! temperature (i.e. zero flux)
120  gw = 1d0
121  hw = 0d0
122  cw = 0d0
123  ELSE
124 ! Setting the wall temperature based on Johnson and Jackson
125 ! Pseudo-thermal temp B.C.
126  CALL calc_theta_bc(ijk,east_of(ijk),'E',&
127  m,l,gw,hw,cw)
128 ! Overwriting calculated value of cw
129  IF (bc_jj_ps(l).EQ.2) cw=0d0
130  ENDIF
131 
132  a_m(ijk, east, m) = -(half*hw - odx_e(i)*gw)
133  a_m(ijk, 0, m) = -(half*hw + odx_e(i)*gw)
134 
135  IF (bc_jj_ps(l) .EQ. 1) THEN
136  b_m(ijk, m) = -cw
137  ELSE
138  b_m(ijk,m) = zero
139  ENDIF
140  ENDIF
141 
142  ELSEIF(fluid_at(west_of(ijk)))THEN
143  IF(ep_s(west_of(ijk),m).LE.dil_ep_s)THEN
144  a_m(ijk, west, m) = one
145  ELSE
146  IF (bc_jj_ps(l).EQ.3) THEN
147  gw = 1d0
148  hw = 0d0
149  cw = 0d0
150  ELSE
151  CALL calc_theta_bc(ijk,west_of(ijk),'W',&
152  m,l,gw,hw,cw)
153  IF (bc_jj_ps(l).EQ.2) cw=0d0
154  ENDIF
155 
156  a_m(ijk, west, m) = -(half*hw - odx_e(im)*gw)
157  a_m(ijk, 0, m) = -(half*hw + odx_e(im)*gw)
158 
159  IF (bc_jj_ps(l) .EQ. 1) THEN
160  b_m(ijk, m) = -cw
161  ELSE
162  b_m(ijk,m) = zero
163  ENDIF
164  ENDIF
165 
166  ELSEIF(fluid_at(north_of(ijk)))THEN
167  IF(ep_s(north_of(ijk),m).LE.dil_ep_s)THEN
168  a_m(ijk, north, m) = one
169  ELSE
170  IF (bc_jj_ps(l).EQ.3) THEN
171  gw = 1d0
172  hw = 0d0
173  cw = 0d0
174  ELSE
175  CALL calc_theta_bc(ijk,north_of(ijk),'N',&
176  m,l,gw,hw,cw)
177  IF (bc_jj_ps(l).EQ.2) cw=0d0
178  ENDIF
179 
180  a_m(ijk, north, m) = -(half*hw - ody_n(j)*gw)
181  a_m(ijk, 0, m) = -(half*hw + ody_n(j)*gw)
182 
183  IF (bc_jj_ps(l) .EQ. 1) THEN
184  b_m(ijk, m) = -cw
185  ELSE
186  b_m(ijk,m) = zero
187  ENDIF
188  ENDIF
189 
190  ELSEIF(fluid_at(south_of(ijk)))THEN
191  IF(ep_s(south_of(ijk),m).LE.dil_ep_s)THEN
192  a_m(ijk, south, m) = one
193  ELSE
194  IF (bc_jj_ps(l).EQ.3) THEN
195  gw = 1d0
196  hw = 0d0
197  cw = 0d0
198  ELSE
199  CALL calc_theta_bc(ijk,south_of(ijk),'S',&
200  m,l,gw,hw,cw)
201  IF (bc_jj_ps(l).EQ.2) cw=0d0
202  ENDIF
203 
204  a_m(ijk, south, m) = -(half*hw - ody_n(jm)*gw)
205  a_m(ijk, 0, m) = -(half*hw + ody_n(jm)*gw)
206 
207  IF (bc_jj_ps(l) .EQ. 1) THEN
208  b_m(ijk, m) = -cw
209  ELSE
210  b_m(ijk,m) = zero
211  ENDIF
212  ENDIF
213 
214  ELSEIF(fluid_at(top_of(ijk)))THEN
215  IF(ep_s(top_of(ijk),m).LE.dil_ep_s)THEN
216  a_m(ijk, top, m) = one
217  ELSE
218  IF (bc_jj_ps(l).EQ.3) THEN
219  gw = 1d0
220  hw = 0d0
221  cw = 0d0
222  ELSE
223  CALL calc_theta_bc(ijk,top_of(ijk),'T',&
224  m,l,gw,hw,cw)
225  IF (bc_jj_ps(l).EQ.2) cw=0d0
226  ENDIF
227 
228  a_m(ijk, top, m) = -(half*hw - ox(i)*odz_t(k)*gw)
229  a_m(ijk, 0, m) = -(half*hw + ox(i)*odz_t(k)*gw)
230 
231  IF (bc_jj_ps(l) .EQ. 1) THEN
232  b_m(ijk, m) = -cw
233  ELSE
234  b_m(ijk,m) = zero
235  ENDIF
236  ENDIF
237 
238  ELSEIF(fluid_at(bottom_of(ijk)))THEN
239  IF(ep_s(bottom_of(ijk),m).LE.dil_ep_s)THEN
240  a_m(ijk, bottom, m) = one
241  ELSE
242  IF (bc_jj_ps(l).EQ.3) THEN
243  gw = 1d0
244  hw = 0d0
245  cw = 0d0
246  ELSE
247  CALL calc_theta_bc(ijk,bottom_of(ijk),'B',&
248  m,l,gw,hw,cw)
249  IF (bc_jj_ps(l).EQ.2) cw=0d0
250  ENDIF
251 
252  a_m(ijk, bottom, m) = -(half*hw - ox(i)*odz_t(km)*gw)
253  a_m(ijk, 0, m) = -(half*hw + ox(i)*odz_t(km)*gw)
254 
255  IF (bc_jj_ps(l) .EQ. 1) THEN
256  b_m(ijk, m) = -cw
257  ELSE
258  b_m(ijk,m) = zero
259  ENDIF
260  ENDIF
261  ENDIF
262 
263  ENDDO ! end do over i
264  ENDDO ! end do over j
265  ENDDO ! end do over k
266 
267  ENDIF ! end if (bc_jj_ps>0)
268  ENDIF ! end if nsw, fsw or psw
269  ENDIF ! end if (bc_defined)
270  ENDDO ! end L do loop over dimension_bc
271 
272  RETURN
273  END SUBROUTINE bc_theta
274 
275 
276 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
277 ! C
278 ! Subroutine: CALC_THETA_BC C
279 ! C
280 ! Purpose: Implementation of Johnson & Jackson boundary conditions C
281 ! for the pseudo-thermal temperature. C
282 ! C
283 ! C
284 ! Author: Kapil Agrawal, Princeton University Date: 14-MAR-98 C
285 ! Reviewer: Date: C
286 ! C
287 ! C
288 ! C
289 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
290 
291  SUBROUTINE calc_theta_bc(IJK1,IJK2,FCELL,M,L,Gw,Hw,Cw)
292 
293 !-----------------------------------------------
294 ! Modules
295 !-----------------------------------------------
296  USE bc
297  USE compar
298  USE constant
299  USE fldvar
300  USE fun_avg
301  USE functions
302  USE geometry
303  USE indices
304  USE mpi_utility
305  USE param
306  USE param1
307  USE physprop
308  USE rdf
309  USE run
310  USE rxns
311  USE toleranc
312  USE turb
313  USE visc_s
314  IMPLICIT NONE
315 !-----------------------------------------------
316 ! Dummy Arguments
317 !-----------------------------------------------
318 ! IJK indices for wall cell (ijk1) and fluid cell (ijk2)
319  INTEGER, INTENT(IN) :: IJK1, IJK2
320 ! The location (e,w,n...) of fluid cell
321  CHARACTER, INTENT(IN) :: FCELL
322 ! Solids phase index
323  INTEGER, INTENT(IN) :: M
324 ! Index corresponding to boundary condition
325  INTEGER, INTENT(IN) :: L
326 ! Granular energy coefficient
327  DOUBLE PRECISION, INTENT(INOUT) :: Gw, Hw, Cw
328 !-----------------------------------------------
329 ! Local Variables
330 !-----------------------------------------------
331 ! IJK indices for cells
332  INTEGER :: IJK3
333 ! Solids phase index
334  INTEGER :: MM
335 ! Average scalars
336  DOUBLE PRECISION :: EPg_avg, Mu_g_avg, RO_g_avg, smallTheta
337 ! void fraction at packing
338  DOUBLE PRECISION :: ep_star_avg
339 ! Average scalars modified to include all solid phases
340  DOUBLE PRECISION :: EPs_avg(dimension_m), DP_avg(dimension_m), &
341  TH_avg(dimension_m)
342 ! Average density of solids phase M
343  DOUBLE PRECISION ROs_avg(dimension_m)
344 ! Average Simonin variables
345  DOUBLE PRECISION :: K_12_avg, Tau_12_avg
346 ! values of U_sm, V_sm, W_sm at appropriate place on boundary wall
347  DOUBLE PRECISION :: USCM, VSCM, WSCM
348 ! values of U_g, V_g, W_g at appropriate place on boundary wall
349  DOUBLE PRECISION :: UGC, VGC, WGC
350 ! Magnitude of gas-solids relative velocity
351  DOUBLE PRECISION :: VREL
352 ! Square of slip velocity between wall and particles
353  DOUBLE PRECISION :: VSLIPSQ
354 ! Wall velocity dot outward normal for all solids phases
355  DOUBLE PRECISION :: VWDOTN (dimension_m)
356 ! Gradient in number density at the wall dot with outward normal for
357 ! all solids phases
358  DOUBLE PRECISION :: GNUWDOTN (dimension_m)
359 ! Gradient in granular temperature at the wall dot with outward normal
360 ! for all solids phases
361  DOUBLE PRECISION :: GTWDOTN (dimension_m)
362 ! Error message
363  CHARACTER(LEN=80) LINE(1)
364 ! Radial distribution function
365  DOUBLE PRECISION :: g0(dimension_m)
366 ! Sum of eps*G_0
367  DOUBLE PRECISION :: g0EPs_avg
368 !-----------------------------------------------
369 ! External functions
370 !-----------------------------------------------
371 ! Variable specularity coefficient
372  DOUBLE PRECISION :: PHIP_JJ
373 !-----------------------------------------------
374 
375 ! Note: EP_s, MU_g, and RO_g are undefined at IJK1 (wall cell).
376 ! Hence IJK2 (fluid cell) is used in averages.
377 
378  smalltheta = (to_si)**4 * zero_ep_s
379 
380 ! set location independent quantities
381  epg_avg = ep_g(ijk2)
382  ep_star_avg = ep_star_array(ijk2)
383  mu_g_avg = mu_g(ijk2)
384  ro_g_avg = ro_g(ijk2)
385  g0eps_avg = zero
386 
387 ! added for Simonin model (sof)
388  IF(kt_type_enum == simonin_1996) THEN
389  k_12_avg = k_12(ijk2)
390  tau_12_avg = tau_12(ijk2)
391  ELSE
392  k_12_avg = zero
393  tau_12_avg = zero
394  ENDIF
395 
396  DO mm = 1, smax
397  g0(mm) = g_0(ijk2, m, mm)
398  eps_avg(mm) = ep_s(ijk2,mm)
399  dp_avg(mm) = d_p(ijk2,mm)
400  g0eps_avg = g0eps_avg + g_0(ijk2, m, mm)*ep_s(ijk2,mm)
401  ros_avg(mm) = ro_s(ijk2,mm)
402  ENDDO
403 
404 
405  IF(fcell .EQ. 'N')THEN
406  DO mm = 1, smax
407  th_avg(mm) = avg_y(theta_m(ijk1,mm),theta_m(ijk2,mm),j_of(ijk1))
408  IF(th_avg(mm) < zero) th_avg(mm) = smalltheta ! for some corner cells
409 
410 ! added for IA (2005) theory:
411 ! include -1 since normal vector is pointing south (-y)
412 ! velocity at wall (i,j+1/2,k relative to ijk1) dot with the normal
413 ! vector at the wall
414  vwdotn(mm) = -1.d0*v_s(ijk1,mm)
415 
416 ! gradient in number density at wall (i,j+1/2,k relative to ijk1) dot
417 ! with the normal vector at the wall. since nu is undefined at ijk1,
418 ! approximate gradient with interior points: ijk2 and i,j+1,k relative
419 ! to ijk2
420  ijk3 = north_of(ijk2)
421  gnuwdotn(mm) = -1.d0*(6.d0/(pi*dp_avg(mm)))*&
422  (ep_s(ijk3,mm)-ep_s(ijk2,mm))*ody_n(j_of(ijk2))
423 
424 ! gradient in granular temperature at wall (i,j+1/2,k relative to ijk1) dot
425 ! with the normal vector at the wall.
426  gtwdotn(mm) = -1.d0*(theta_m(ijk3,mm)-theta_m(ijk2,mm))*&
427  ody_n(j_of(ijk2))
428  ENDDO
429 
430 ! Calculate velocity components at i, j+1/2, k (relative to IJK1)
431  ugc = avg_y(avg_x_e(u_g(im_of(ijk1)),u_g(ijk1),i_of(ijk1)),&
432  avg_x_e(u_g(im_of(ijk2)),u_g(ijk2),i_of(ijk2)),&
433  j_of(ijk1))
434  vgc = v_g(ijk1)
435  wgc = avg_y(avg_z_t(w_g(km_of(ijk1)), w_g(ijk1)),&
436  avg_z_t(w_g(km_of(ijk2)), w_g(ijk2)),&
437  j_of(ijk1))
438  uscm = avg_y(avg_x_e(u_s(im_of(ijk1),m),u_s(ijk1,m),i_of(ijk1)),&
439  avg_x_e(u_s(im_of(ijk2),m),u_s(ijk2,m),i_of(ijk2)),&
440  j_of(ijk1))
441  vscm = v_s(ijk1,m)
442  wscm = avg_y(avg_z_t(w_s(km_of(ijk1),m), w_s(ijk1,m)),&
443  avg_z_t(w_s(km_of(ijk2),m), w_s(ijk2,m)),&
444  j_of(ijk1))
445 
446 ! slip velocity at the wall
447  vslipsq=(wscm-bc_ww_s(l, m))**2 + (uscm-bc_uw_s(l, m))**2
448 
449 
450  ELSEIF(fcell .EQ. 'S')THEN
451  DO mm = 1, smax
452  th_avg(mm) = avg_y(theta_m(ijk2, mm),theta_m(ijk1, mm),j_of(ijk2))
453  IF(th_avg(mm) < zero) th_avg(mm) = smalltheta ! for some corner cells
454 
455 ! added for IA (2005) theory:
456 ! include 1 since normal vector is pointing north (+y)
457 ! velocity at wall (i,j+1/2,k relative to ijk2) dot with the normal
458 ! vector at the wall.
459  vwdotn(mm) = 1.d0*v_s(ijk2,mm)
460 
461 ! gradient in number density at wall (i,j+1/2,k relative to ijk2) dot
462 ! with the normal vector at the wall. since nu is undefined at ijk1,
463 ! approximate gradient with interior points: ijk2 and i,j-1,k relative
464 ! to ijk2
465  ijk3 = south_of(ijk2)
466  gnuwdotn(mm) = 1.d0*(6.d0/(pi*dp_avg(mm)))*&
467  (ep_s(ijk2,mm)-ep_s(ijk3,mm))*ody_n(j_of(ijk3))
468 
469 ! gradient in granular temperature at wall (i,j+1/2,k relative to ijk2)
470 ! dot with the normal vector at the wall.
471  gtwdotn(mm) = 1.d0*(theta_m(ijk2,mm)-theta_m(ijk3,mm))*&
472  ody_n(j_of(ijk3))
473  ENDDO
474 
475 ! Calculate velocity components at i, j+1/2, k (relative to IJK2)
476  ugc = avg_y(avg_x_e(u_g(im_of(ijk2)),u_g(ijk2),i_of(ijk2)),&
477  avg_x_e(u_g(im_of(ijk1)),u_g(ijk1),i_of(ijk1)),&
478  j_of(ijk2))
479  vgc = v_g(ijk2)
480  wgc = avg_y(avg_z_t(w_g(km_of(ijk2)), w_g(ijk2)),&
481  avg_z_t(w_g(km_of(ijk1)), w_g(ijk1)),&
482  j_of(ijk2))
483  uscm = avg_y(avg_x_e(u_s(im_of(ijk2),m),u_s(ijk2,m),i_of(ijk2)),&
484  avg_x_e(u_s(im_of(ijk1),m),u_s(ijk1,m),i_of(ijk1)),&
485  j_of(ijk2))
486  vscm = v_s(ijk2,m)
487  wscm = avg_y(avg_z_t(w_s(km_of(ijk2),m), w_s(ijk2,m)),&
488  avg_z_t(w_s(km_of(ijk1),m), w_s(ijk1,m)),&
489  j_of(ijk2))
490 
491 ! slip velocity at the wall
492  vslipsq=(wscm-bc_ww_s(l, m))**2 + (uscm-bc_uw_s(l, m))**2
493 
494 
495  ELSEIF(fcell== 'E')THEN
496  DO mm = 1, smax
497  th_avg(mm) = avg_x(theta_m(ijk1,mm),theta_m(ijk2,mm),i_of(ijk1))
498  IF(th_avg(mm) < zero) th_avg(mm) = smalltheta ! for some corner cells
499 
500 ! added for IA (2005) theory:
501 ! include -1 since normal vector is pointing west (-x)
502 ! velocity at wall (i+1/2,j,k relative to ijk1) dot with the normal
503 ! vector at the wall.
504  vwdotn(mm) = -1.d0*u_s(ijk1,mm)
505 
506 ! gradient in number density at wall (i+1/2,j,k relative to ijk1) dot
507 ! with the normal vector at the wall. since nu is undefined at ijk1,
508 ! approximate gradient with interior points: ijk2 and i,i+1,k relative
509 ! to ijk2
510  ijk3 = east_of(ijk2)
511  gnuwdotn(mm) = -1.d0*(6.d0/(pi*dp_avg(mm)))*&
512  (ep_s(ijk3,mm)-ep_s(ijk2,mm))*odx_e(i_of(ijk2))
513 
514 ! gradient in granular temperature at wall (i+1/2,j,k relative to ijk1)
515 ! dot with the normal vector at the wall.
516  gtwdotn(mm) = -1.d0*(theta_m(ijk3,mm)-theta_m(ijk2,mm))*&
517  odx_e(i_of(ijk2))
518  ENDDO
519 
520 ! Calculate velocity components at i+1/2, j, k (relative to IJK1)
521  ugc = u_g(ijk1)
522  vgc = avg_x(avg_y_n(v_g(jm_of(ijk1)), v_g(ijk1)),&
523  avg_y_n(v_g(jm_of(ijk2)), v_g(ijk2)),&
524  i_of(ijk1))
525  wgc = avg_x(avg_z_t(w_g(km_of(ijk1)), w_g(ijk1)),&
526  avg_z_t(w_g(km_of(ijk2)), w_g(ijk2)),&
527  i_of(ijk1))
528  uscm = u_s(ijk1,m)
529  vscm = avg_x(avg_y_n(v_s(jm_of(ijk1),m), v_s(ijk1,m)),&
530  avg_y_n(v_s(jm_of(ijk2),m), v_s(ijk2,m)),&
531  i_of(ijk1))
532  wscm = avg_x(avg_z_t(w_s(km_of(ijk1),m), w_s(ijk1,m)),&
533  avg_z_t(w_s(km_of(ijk2),m), w_s(ijk2,m)),&
534  i_of(ijk1))
535 
536 ! slip velocity at the wall
537  vslipsq=(wscm-bc_ww_s(l, m))**2 + (vscm-bc_vw_s(l, m))**2
538 
539 
540  ELSEIF(fcell== 'W')THEN
541  DO mm = 1, smax
542  th_avg(mm) = avg_x(theta_m(ijk2,mm),theta_m(ijk1,mm),i_of(ijk2))
543  IF(th_avg(mm) < zero) th_avg(mm) = smalltheta ! for some corner cells
544 
545 ! added for IA (2005) theory:
546 ! include 1 since normal vector is pointing west (+x)
547 ! velocity at wall (i+1/2,j,k relative to ijk2) dot with the normal
548 ! vector at the wall.
549  vwdotn(mm) = 1.d0*u_s(ijk2,mm)
550 
551 ! gradient in number density at wall (i+1/2,j,k relative to ijk2) dot
552 ! with the normal vector at the wall. since nu is undefined at ijk1,
553 ! approximate gradient with interior points: ijk2 and i,i-1,k relative
554 ! to ijk2
555  ijk3 = west_of(ijk2)
556  gnuwdotn(mm) = 1.d0*(6.d0/(pi*dp_avg(mm)))*&
557  (ep_s(ijk2,mm)-ep_s(ijk3,mm))*odx_e(i_of(ijk3))
558 
559 ! gradient in granular temperature at wall (i+1/2,j,k relative to ijk2)
560 ! dot with the normal vector at the wall.
561  gtwdotn(mm) = 1.d0*(theta_m(ijk2,mm)-theta_m(ijk3,mm))*&
562  odx_e(i_of(ijk3))
563  ENDDO
564 
565 ! Calculate velocity components at i+1/2, j, k (relative to IJK2)
566  ugc = u_g(ijk2)
567  vgc = avg_x(avg_y_n(v_g(jm_of(ijk2)), v_g(ijk2)),&
568  avg_y_n(v_g(jm_of(ijk1)), v_g(ijk1)),&
569  i_of(ijk2))
570  wgc = avg_x(avg_z_t(w_g(km_of(ijk2)), w_g(ijk2)),&
571  avg_z_t(w_g(km_of(ijk1)), w_g(ijk1)),&
572  i_of(ijk2))
573  uscm = u_s(ijk2,m)
574  vscm = avg_x(avg_y_n(v_s(jm_of(ijk2),m), v_s(ijk2,m)),&
575  avg_y_n(v_s(jm_of(ijk1),m), v_s(ijk1,m)),&
576  i_of(ijk2))
577  wscm = avg_x(avg_z_t(w_s(km_of(ijk2),m), w_s(ijk2,m)),&
578  avg_z_t(w_s(km_of(ijk1),m), w_s(ijk1,m)),&
579  i_of(ijk2))
580 
581 ! slip velocity at the wall
582  vslipsq=(wscm-bc_ww_s(l, m))**2 + (vscm-bc_vw_s(l, m))**2
583 
584 
585  ELSEIF(fcell== 'T')THEN
586  DO mm = 1, smax
587  th_avg(mm) = avg_z(theta_m(ijk1,mm),theta_m(ijk2,mm),k_of(ijk1))
588  IF(th_avg(mm) < zero) th_avg(mm) = smalltheta ! for some corner cells
589 
590 ! added for IA (2005) theory:
591 ! include -1 since normal vector is pointing in bottom dir (-z)
592 ! velocity at wall (i,j,k+1/2 relative to ijk1) dot with the normal
593 ! vector at the wall.
594  vwdotn(mm) = -1.d0*w_s(ijk1,mm)
595 
596 ! gradient in number density at wall (i,j,k+1/2 relative to ijk1) dot
597 ! with the normal vector at the wall. since nu is undefined at ijk1,
598 ! approximate gradient with interior points: ijk2 and i,j,k+1 relative
599 ! to ijk2
600  ijk3 = top_of(ijk2)
601  gnuwdotn(mm) = -1.d0*(6.d0/(pi*dp_avg(mm)))*&
602  (ep_s(ijk3,mm)-ep_s(ijk2,mm))*ox(i_of(ijk2))*odz_t(k_of(ijk2))
603 
604 ! gradient in granular temperature at wall (i,j,k+1/2 relative to ijk1)
605 ! dot with the normal vector at the wall.
606  gnuwdotn(mm) = -1.d0*(theta_m(ijk3,mm)-theta_m(ijk2,mm))*&
607  ox(i_of(ijk2))*odz_t(k_of(ijk2))
608  ENDDO
609 
610 ! Calculate velocity components at i, j, k+1/2 (relative to IJK1)
611  ugc = avg_z(avg_x_e(u_g(im_of(ijk1)),u_g(ijk1),i_of(ijk1)),&
612  avg_x_e(u_g(im_of(ijk2)),u_g(ijk2),i_of(ijk2)),&
613  k_of(ijk1))
614  vgc = avg_z(avg_y_n(v_g(jm_of(ijk1)), v_g(ijk1)),&
615  avg_y_n(v_g(jm_of(ijk2)), v_g(ijk2)),&
616  k_of(ijk1))
617  wgc = w_g(ijk1)
618  uscm = avg_z(avg_x_e(u_s(im_of(ijk1),m),u_s(ijk1,m),i_of(ijk1)),&
619  avg_x_e(u_s(im_of(ijk2),m),u_s(ijk2,m),i_of(ijk2)),&
620  k_of(ijk1))
621  vscm = avg_z(avg_y_n(v_s(jm_of(ijk1),m), v_s(ijk1,m)),&
622  avg_y_n(v_s(jm_of(ijk2),m), v_s(ijk2,m)),&
623  k_of(ijk1))
624  wscm = w_s(ijk1,m)
625 
626 ! slip velocity at the wall
627  vslipsq=(vscm-bc_vw_s(l, m))**2 + (uscm-bc_uw_s(l, m))**2
628 
629 
630  ELSEIF(fcell== 'B')THEN
631  DO mm = 1, smax
632  th_avg(mm) = avg_z(theta_m(ijk2,mm),theta_m(ijk1,mm),k_of(ijk2))
633  IF(th_avg(mm) < zero) th_avg(mm) = smalltheta ! for some corner cells
634 
635 ! added for IA (2005) theory:
636 ! include 1 since normal vector is pointing in bottom dir (+z)
637 ! velocity at wall (i,j,k+1/2 relative to ijk2) dot with the normal
638 ! vector at the wall
639  vwdotn(mm) = 1.d0*w_s(ijk2,mm)
640 
641 ! gradient in number density at wall (i,j,k+1/2 relative to ijk2) dot
642 ! with the normal vector at the wall. since nu is undefined at ijk1,
643 ! approximate gradient with interior points: ijk2 and i,j,k-1 relative
644 ! to ijk2
645  ijk3 = bottom_of(ijk2)
646  gnuwdotn(mm) = 1.d0*(6.d0/(pi*dp_avg(mm)))*&
647  (ep_s(ijk2,mm)-ep_s(ijk3,mm))*ox(i_of(ijk3))*odz_t(k_of(ijk3))
648 
649 ! gradient in granular temperature at wall (i,j,k+1/2 relative to ijk2)
650 ! dot with the normal vector at the wall.
651  gtwdotn(mm) = 1.d0*(theta_m(ijk2,mm)-theta_m(ijk3,mm))*&
652  ox(i_of(ijk3))*odz_t(k_of(ijk3))
653  ENDDO
654 
655 ! Calculate velocity components at i, j, k+1/2 (relative to IJK2)
656  ugc = avg_z(avg_x_e(u_g(im_of(ijk2)),u_g(ijk2),i_of(ijk2)),&
657  avg_x_e(u_g(im_of(ijk1)),u_g(ijk1),i_of(ijk1)),&
658  k_of(ijk2))
659  vgc = avg_z(avg_y_n(v_g(jm_of(ijk2)), v_g(ijk2)),&
660  avg_y_n(v_g(jm_of(ijk1)), v_g(ijk1)),&
661  k_of(ijk2))
662  wgc = w_g(ijk2)
663  uscm = avg_z(avg_x_e(u_s(im_of(ijk2),m),u_s(ijk2,m),i_of(ijk2)),&
664  avg_x_e(u_s(im_of(ijk1),m),u_s(ijk1,m),i_of(ijk1)),&
665  k_of(ijk2))
666  vscm = avg_z(avg_y_n(v_s(jm_of(ijk2),m), v_s(ijk2,m)),&
667  avg_y_n(v_s(jm_of(ijk1),m), v_s(ijk1,m)),&
668  k_of(ijk2))
669  wscm = w_s(ijk2,m)
670 
671 ! slip velocity at the wall
672  vslipsq=(vscm-bc_vw_s(l, m))**2 + (uscm-bc_uw_s(l, m))**2
673 
674  ELSE
675  WRITE(line,'(A, A)') 'Error: Unknown FCELL'
676  CALL write_error('CALC_THETA_BC', line, 1)
677  call exitmpi(mype)
678  ENDIF
679 
680 ! magnitude of gas-solids relative velocity
681  vrel = dsqrt( (ugc-uscm)**2 + (vgc-vscm)**2 + (wgc-wscm)**2 )
682 
683  CALL theta_hw_cw(g0, eps_avg, epg_avg, ep_star_avg, &
684  g0eps_avg, th_avg,mu_g_avg,ro_g_avg, ros_avg, &
685  dp_avg, k_12_avg,tau_12_avg,vrel,vslipsq,vwdotn,&
686  gnuwdotn,gtwdotn,m,gw,hw,cw,l)
687 
688 ! Output the variable specularity coefficient. Currently only works
689 ! for one solids phase
690  IF(bc_jj_m .and. phip_out_jj) THEN
691  IF(phip_out_iter.eq.1) THEN
692  reactionrates(ijk1,1)= phip_jj(dsqrt(vslipsq),th_avg(1))
693  ENDIF
694  ENDIF
695 
696  RETURN
697  END SUBROUTINE calc_theta_bc
698 
699 
700 
701 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
702 ! C
703 ! Subroutine: SUBROUTINE THETA_HW_CW C
704 ! Purpose: Subroutine for gw, hw and cw C
705 ! C
706 ! Author: Kapil Agrawal, Princeton University Date: 15-MAR-98 C
707 ! Reviewer: Date: C
708 ! C
709 ! C
710 ! Modified: Sofiane Benyahia, Fluent Inc. Date: 02-FEB-05 C
711 ! Purpose: Include conductivity defined by Simonin and Ahmadi C
712 ! Also included Jenkins small frictional limit C
713 ! C
714 ! Literature/Document References: C
715 ! See calc_mu_s.f for references on kinetic theory models C
716 ! Johnson, P. C., and Jackson, R., Frictional-collisional C
717 ! constitutive relations for granular materials, with C
718 ! application to plane shearing, Journal of Fluid Mechanics, C
719 ! 1987, 176, pp. 67-93 C
720 ! Jenkins, J. T., and Louge, M. Y., On the flux of fluctuating C
721 ! energy in a collisional grain flow at a flat frictional wall, C
722 ! Physics of Fluids, 1997, 9(10), pp. 2835-2840 C
723 ! C
724 ! C
725 ! Additional Notes: C
726 ! The JJ granular energy BC is written as: n.q = D-G, where C
727 ! q = the heat flux, n = outward normal vector C
728 ! D = dissipation due to p-w collisions, and C
729 ! G = generation due to p-w slip C
730 ! In MFIX this is implemented as k.gradT = G-D, where C
731 ! k = granular conductivity units (Mass/Length/Time) C
732 ! T = granular temperature in units (Length/Time)^2 C
733 ! However, the expression for heat flux (q) may contain terms C
734 ! beyond that of gradient in temperature, such as the dufour C
735 ! term. A more rigorous implemenation should account for these C
736 ! additional terms, which is not done here. C
737 ! C
738 ! C
739 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
740 
741  SUBROUTINE theta_hw_cw(g0,EPS,EPG, ep_star_avg, &
742  g0eps_avg,th,mu_g_avg,ro_g_avg, ros_avg, dp_avg, &
743  k_12_avg,tau_12_avg,vrel,vslipsq,vwdotn,&
744  gnuwdotn,gtwdotn,m,gw,hw,cw,l)
745 
746 !-----------------------------------------------
747 ! Modules
748 !-----------------------------------------------
749  USE bc
750  USE compar
751  USE constant
752  USE exit, only: mfix_exit
753  USE fldvar
754  USE kintheory
755  USE param
756  USE param1
757  USE physprop
758  USE run
759  USE toleranc
760  IMPLICIT NONE
761 !-----------------------------------------------
762 ! Dummy Arguments
763 !-----------------------------------------------
764 ! Radial distribution function of solids phase M with each
765 ! other solids phase
766  DOUBLE PRECISION, INTENT(IN) :: g0(dimension_m)
767 ! Average solids volume fraction of each solids phase
768  DOUBLE PRECISION, INTENT(IN) :: EPS(dimension_m)
769 ! Average solids and gas volume fraction
770  DOUBLE PRECISION, INTENT(IN) :: EPG, ep_star_avg
771 ! Sum of eps*G_0
772  DOUBLE PRECISION, INTENT(INOUT) :: g0EPs_avg
773 ! Average theta_m
774  DOUBLE PRECISION, INTENT(INOUT) :: TH (dimension_m)
775 ! Average gas viscosity
776  DOUBLE PRECISION, INTENT(IN) :: Mu_g_avg
777 ! Average gas density
778  DOUBLE PRECISION, INTENT(IN) :: RO_g_avg
779 ! Average density of solids phase M
780  DOUBLE PRECISION, INTENT(IN) :: ROS_avg(dimension_m)
781 ! Average particle diameter of each solids phase
782  DOUBLE PRECISION, INTENT(IN) :: DP_avg(dimension_m)
783 ! Average cross-correlation K_12 and lagrangian integral time-scale
784  DOUBLE PRECISION, INTENT(IN) :: K_12_avg, Tau_12_avg
785 ! Magnitude of slip velocity between two phases
786  DOUBLE PRECISION, INTENT(IN) :: VREL
787 ! Square of slip velocity between wall and particles
788  DOUBLE PRECISION, INTENT(IN) :: VSLIPSQ
789 ! Wall velocity dot outward normal for all solids phases
790  DOUBLE PRECISION, INTENT(IN) :: VWDOTN (dimension_m)
791 ! Gradient in number density at the wall dot with outward
792 ! normal for all solids phases
793  DOUBLE PRECISION, INTENT(IN) :: GNUWDOTN (dimension_m)
794 ! Gradient in temperature at the wall dot with outward
795 ! normal for all solids phases
796  DOUBLE PRECISION, INTENT(IN) :: GTWDOTN (dimension_m)
797 ! Solids phase index
798  INTEGER, INTENT(IN) :: M
799 ! Coefficients in JJ boundary conditions
800 ! Gw = granular conductivity (w particle mass)
801 ! Hw = dissipation due to particle-wall collision
802 ! Cw = generation due to particle-wall slip
803  DOUBLE PRECISION, INTENT(INOUT) :: GW, HW, CW
804 ! Index corresponding to boundary condition
805  INTEGER, INTENT(IN) :: L
806 !-----------------------------------------------
807 ! Local Variables
808 !-----------------------------------------------
809 ! Solids phase index
810  INTEGER :: LL
811 ! Coefficient of 1st term
812  DOUBLE PRECISION :: K_1
813 ! Conductivity
814  DOUBLE PRECISION :: Kgran
815 ! Conductivity corrected for interstitial fluid effects
816  DOUBLE PRECISION :: Kgran_star
817 ! Drag Coefficient
818  DOUBLE PRECISION :: Beta, DgA
819 ! Reynolds number based on slip velocity
820  DOUBLE PRECISION :: Re_g
821 ! Friction Factor in drag coefficient
822  DOUBLE PRECISION :: C_d
823 ! mth solids phase particle diameter, mass and number density
824  DOUBLE PRECISION :: D_PM, M_PM, NU_PM
825 ! Variables for Iddir & Arastoopour model
826  DOUBLE PRECISION :: D_PL, M_PL, NU_PL, MPSUM, DPSUMo2
827  DOUBLE PRECISION :: Ap_lm, Dp_lm, R0p_lm, R1p_lm, R8p_lm, R9p_lm, &
828  Bp_lm, R5p_lm, R6p_lm, R7p_lm
829  DOUBLE PRECISION :: K_s_sum, K_s_MM, K_sM_LM
830  DOUBLE PRECISION :: Kvel_s_sum, Kth_s_sum
831  DOUBLE PRECISION :: Knu_s_sum, K_common_term
832 ! Variables for Garzo & Dufty model
833  DOUBLE PRECISION :: c_star, zeta0_star, press_star, eta0, &
834  kappa0, nu_kappa_star, kappa_k_star, &
835  kappa_star
836 ! Variables for GTSH theory
837  DOUBLE PRECISION :: Re_T, Chi, vfrac, mu2_0, mu4_0, mu4_1, &
838  A2_gtshW, zeta_star, nu0, NuK, Kth0, KthK, EDT_s
839 ! Local parameters for Simonin model
840  DOUBLE PRECISION :: Zeta_c, Omega_c, Tau_2_c, Kappa_kin, &
841  Kappa_Col, Tau_12_st
842 ! Solids pressure and pressure divided by granular temperature
843  DOUBLE PRECISION :: Ps, PsoTheta
844 !-----------------------------------------------
845 ! Function subroutines
846 !-----------------------------------------------
847 ! variable specularity coefficient
848  DOUBLE PRECISION :: PHIP_JJ
849 !-----------------------------------------------
850 
851  IF(th(m) .LE. zero)THEN
852  th(m) = 1d-8
853  IF (mype.eq.pe_io) THEN
854  WRITE(*,*) &
855  'Warning: Negative granular temp at wall set to 1e-8'
856 ! CALL WRITE_ERROR('THETA_HW_CW', LINE, 1)
857  ENDIF
858  ENDIF
859 
860 ! common variables
861  d_pm = dp_avg(m)
862  m_pm = (pi/6.d0)*(d_pm**3.)*ros_avg(m)
863  nu_pm = (eps(m)*ros_avg(m))/m_pm
864 
865 ! This is from Wen-Yu correlation, you can put here your own single particle drag
866  re_g = epg*ro_g_avg*d_pm*vrel/mu_g_avg
867  IF (re_g .lt. 1000.d0) THEN
868  c_d = (24.d0/(re_g+small_number))*(one + 0.15d0 * re_g**0.687d0)
869  ELSE
870  c_d = 0.44d0
871  ENDIF
872  dga = 0.75d0*c_d*ro_g_avg*epg*vrel/(dp_avg(m)*epg**(2.65d0))
873  IF(vrel == zero) dga = large_number
874  beta = eps(m)*dga
875 
876 
877 ! Defining conductivity according to each KT for later use in JJ BC
878 ! Also define solids pressure for use in Jenkins small frictional BC
879 ! -------------------------------------------------------------------
880  SELECT CASE (kt_type_enum)
881 
882  CASE (lun_1984)
883  kgran = 75d0*ros_avg(m)*dp_avg(m)*dsqrt(pi*th(m))/&
884  (48*eta*(41d0-33d0*eta))
885 
886  IF(switch == zero .OR. ro_g_avg == zero)THEN
887  kgran_star = kgran
888  ELSEIF(th(m) .LT. small_number)THEN
889  kgran_star = zero
890  ELSE
891  kgran_star = ros_avg(m)*eps(m)* g0(m)*th(m)* kgran/ &
892  (ros_avg(m)*g0eps_avg*th(m) + &
893  1.2d0*dga/ros_avg(m)* kgran)
894  ENDIF
895 
896 ! mth solids phase granular conductivity includes particle mass
897  k_1 = kgran_star/g0(m)*( &
898  ( one + (12d0/5.d0)*eta*g0eps_avg ) * &
899  ( one + (12d0/5.d0)*eta*eta*(4d0*eta-3d0) *g0eps_avg ) + &
900  (64d0/(25d0*pi)) * (41d0-33d0*eta) *(eta*g0eps_avg)**2 )
901 
902 ! defining granular pressure (for Jenkins BC)
903  psotheta = ros_avg(m)*eps(m)*(one+4.d0*eta*g0eps_avg)
904  ps = ros_avg(m)*eps(m)*th(m)*(one+4.d0*eta*g0eps_avg)
905 
906 
907  CASE (simonin_1996)
908 ! particle relaxation time
909  tau_12_st = ros_avg(m)/(dga+small_number)
910  zeta_c = (one+ c_e)*(49.d0-33.d0*c_e)/100.d0
911  omega_c = 3.d0*(one+ c_e)**2 *(2.d0*c_e-one)/5.d0
912  tau_2_c = dp_avg(m)/(6.d0*eps(m)*g0(m) &
913  *dsqrt(16.d0*(th(m)+small_number)/pi))
914 
915 ! Defining Simonin's Solids Turbulent Kinetic diffusivity: Kappa
916  kappa_kin = (9.d0/10.d0*k_12_avg*(tau_12_avg/tau_12_st) &
917  + 3.0d0/2.0d0 * th(m)*(one+ omega_c*eps(m)*g0(m)))/ &
918  (9.d0/(5.d0*tau_12_st) + zeta_c/tau_2_c)
919 
920  kappa_col = 18.d0/5.d0*eps(m)*g0(m)*eta* (kappa_kin+ &
921  5.d0/9.d0*dp_avg(m)*dsqrt(th(m)/pi))
922 
923 ! mth solids phase granular conductivity includes particle mass
924  k_1 = eps(m)*ros_avg(m)*(kappa_kin + kappa_col)
925 
926 ! defining granular pressure (for Jenkins BC)
927  psotheta = ros_avg(m)*eps(m)*(one+4.d0*eta*g0eps_avg)
928  ps = ros_avg(m)*eps(m)*th(m)*(one+4.d0*eta*g0eps_avg)
929 
930 
931  CASE (ahmadi_1995)
932 ! mth solids phase granular conductivity includes particle mass
933  k_1 = 0.1306d0*ros_avg(m)*dp_avg(m)*(one+c_e**2)* ( &
934  one/g0(m)+4.8d0*eps(m)+12.1184d0 *eps(m)*eps(m)*g0(m) )*&
935  dsqrt(th(m))
936 
937 ! defining granular pressure (for Jenkins BC)
938  psotheta = ros_avg(m)*eps(m)* &
939  ((one + 4.0d0*g0eps_avg)+half*(one -c_e*c_e))
940  ps = ros_avg(m)*eps(m)*th(m)*&
941  ((one + 4.0d0*g0eps_avg)+half*(one -c_e*c_e))
942 
943 
944  CASE (ia_2005)
945 ! Use original IA theory if SWITCH_IA is false
946  IF(.NOT. switch_ia) g0eps_avg = eps(m)*ros_avg(m)
947 
948  k_s_sum = zero
949  kvel_s_sum = zero
950  knu_s_sum = zero
951  kth_s_sum = zero
952 
953  kgran = (75.d0/384.d0)*ros_avg(m)*d_pm*dsqrt(pi*th(m)/m_pm)
954 
955  IF(switch == zero .OR. ro_g_avg == zero)THEN
956  kgran_star = kgran
957  ELSEIF(th(m)/m_pm .LT. small_number)THEN
958  kgran_star = zero
959  ELSE
960  kgran_star = kgran*g0(m)*eps(m)/ &
961  (g0eps_avg+ 1.2d0*dga*kgran / (ros_avg(m)**2 *(th(m)/m_pm)))
962  ENDIF
963 
964  k_s_mm = (kgran_star/(m_pm*g0(m)))*& ! Kth doesn't include the mass.
965  (1.d0+(3.d0/5.d0)*(1.d0+c_e)*(1.d0+c_e)*g0eps_avg)**2
966 
967  DO ll = 1, smax
968  d_pl = dp_avg(ll)
969  m_pl = (pi/6.d0)*(d_pl**3.)*ros_avg(ll)
970  mpsum = m_pm + m_pl
971  dpsumo2 = (d_pm+d_pl)/2.d0
972  nu_pl = (eps(ll)*ros_avg(ll))/m_pl
973 
974  IF ( ll .eq. m) THEN
975  k_s_sum = k_s_sum + k_s_mm
976 ! Kth_sL_LM is zero when LL=M because it cancels with terms from K_s_LM
977 ! Kvel_s_LM should be zero when LL=M (same solids phase)
978 ! Knu_s_LM should be zero when LL=M (same solids phase)
979  ELSE
980  ap_lm = (m_pm*th(ll)+m_pl*th(m))/(2.d0)
981  bp_lm = (m_pm*m_pl*(th(ll)-th(m) ))/(2.d0*mpsum)
982  dp_lm = (m_pl*m_pm*(m_pm*th(m)+m_pl*th(ll) ))/&
983  (2.d0*mpsum*mpsum)
984  r0p_lm = ( 1.d0/( ap_lm**1.5 * dp_lm**2.5 ) )+ &
985  ( (15.d0*bp_lm*bp_lm)/( 2.d0* ap_lm**2.5 * dp_lm**3.5 ) )+&
986  ( (175.d0*(bp_lm**4))/( 8.d0*ap_lm**3.5 * dp_lm**4.5 ) )
987  r1p_lm = ( 1.d0/( (ap_lm**1.5)*(dp_lm**3) ) )+ &
988  ( (9.d0*bp_lm*bp_lm)/( ap_lm**2.5 * dp_lm**4 ) )+&
989  ( (30.d0*bp_lm**4) /( 2.d0*ap_lm**3.5 * dp_lm**5 ) )
990  r5p_lm = ( 1.d0/( ap_lm**2.5 * dp_lm**3 ) )+ &
991  ( (5.d0*bp_lm*bp_lm)/( ap_lm**3.5 * dp_lm**4 ) )+&
992  ( (14.d0*bp_lm**4)/( ap_lm**4.5 * dp_lm**5 ) )
993  r6p_lm = ( 1.d0/( ap_lm**3.5 * dp_lm**3 ) )+ &
994  ( (7.d0*bp_lm*bp_lm)/( ap_lm**4.5 * dp_lm**4 ) )+&
995  ( (126.d0*bp_lm**4)/( 5.d0*ap_lm**5.5 * dp_lm**5 ) )
996  r7p_lm = ( 3.d0/( 2.d0*ap_lm**2.5 * dp_lm**4 ) )+ &
997  ( (10.d0*bp_lm*bp_lm)/( ap_lm**3.5 * dp_lm**5 ) )+&
998  ( (35.d0*bp_lm**4)/( ap_lm**4.5 * dp_lm**6 ) )
999  r8p_lm = ( 1.d0/( 2.d0*ap_lm**1.5 * dp_lm**4 ) )+ &
1000  ( (6.d0*bp_lm*bp_lm)/( ap_lm**2.5 * dp_lm**5 ) )+&
1001  ( (25.d0*bp_lm**4)/( ap_lm**3.5 * dp_lm**6 ) )
1002  r9p_lm = ( 1.d0/( ap_lm**2.5 * dp_lm**3 ) )+ &
1003  ( (15.d0*bp_lm*bp_lm)/( ap_lm**3.5 * dp_lm**4 ) )+&
1004  ( (70.d0*bp_lm**4)/( ap_lm**4.5 * dp_lm**5 ) )
1005  k_common_term = dpsumo2**3 * m_pl*m_pm/(2.d0*mpsum)*&
1006  (1.d0+c_e)*g0(ll) * (m_pm*m_pl)**1.5
1007 
1008 ! solids phase 'conductivity' associated with the
1009 ! gradient in granular temperature of species M
1010  k_sm_lm = - k_common_term*nu_pm*nu_pl*(&
1011  ((dpsumo2*dsqrt(pi)/16.d0)*(3.d0/2.d0)*bp_lm*r5p_lm)+&
1012  ((m_pl/(8.d0*mpsum))*(1.d0-c_e)*(dpsumo2*pi/6.d0)*&
1013  (3.d0/2.d0)*r1p_lm)-(&
1014  ((dpsumo2*dsqrt(pi)/16.d0)*(m_pm/8.d0)*bp_lm*r6p_lm)+&
1015  ((m_pl/(8.d0*mpsum))*(1.d0-c_e)*(dpsumo2*dsqrt(pi)/&
1016  8.d0)*m_pm*r9p_lm)+&
1017  ((dpsumo2*dsqrt(pi)/16.d0)*(m_pl*m_pm/(mpsum*mpsum))*&
1018  m_pl*bp_lm*r7p_lm)+&
1019  ((m_pl/(8.d0*mpsum))*(1.d0-c_e)*(dpsumo2*dsqrt(pi)/&
1020  2.d0)*(m_pm/mpsum)**2 * m_pl*r8p_lm)+&
1021  ((dpsumo2*dsqrt(pi)/16.d0)*(m_pm*m_pl/(2.d0*mpsum))*&
1022  r9p_lm)-&
1023  ((m_pl/(8.d0*mpsum))*(1.d0-c_e)*dpsumo2*dsqrt(pi)*&
1024  (m_pm*m_pl/mpsum)*bp_lm*r7p_lm) )*th(ll) )*&
1025  (th(m)**2 * th(ll)**3)
1026 
1027 ! These lines were commented because they are not currently used
1028 ! solids phase 'conductivity' associated with the gradient in granular
1029 ! temperature of species L
1030 ! Kth_sL_LM = K_common_term*NU_PM*NU_PL*(&
1031 ! (-(DPSUMo2*DSQRT(PI)/16.d0)*(3.d0/2.d0)*Bp_lm*R5p_lm)+&
1032 ! (-(M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*PI/6.d0)*&
1033 ! (3.d0/2.d0)*R1p_lm)+(&
1034 ! ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PL/8.d0)*Bp_lm*R6p_lm)+&
1035 ! ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1036 ! 8.d0)*M_PL*R9p_lm)+&
1037 ! ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PL*M_PM/(MPSUM*MPSUM))*&
1038 ! M_PM*Bp_lm*R7p_lm)+&
1039 ! ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1040 ! 2.d0)*(M_PM*M_PM/(MPSUM*MPSUM))*M_PM*R8p_lm)+&
1041 ! ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PM*M_PL/(2.d0*MPSUM))*&
1042 ! R9p_lm)+&
1043 ! ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*DPSUMo2*DSQRT(PI)*&
1044 ! (M_PM*M_PL/MPSUM)*Bp_lm*R7p_lm) )*TH(M) )*&
1045 ! (TH(LL)*TH(LL)*(TH(M)**3.))*(GTWDOTN(LL))
1046 
1047 ! solids phase 'conductivity' associated with the difference in velocities
1048 ! Kvel_s_LM = K_common_term*NU_PM*NU_PL*&
1049 ! (M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(3.d0*PI/10.d0)*&
1050 ! R0p_lm*( (TH(M)*TH(LL))**(5./2.) )*&
1051 ! (VWDOTN(M)-VWDOTN(LL))
1052 
1053 ! solids phase 'conductivity' associated with the difference in the
1054 ! gradient in number densities
1055 ! Knu_s_LM = K_common_term*(((DPSUMo2*DSQRT(PI)/16.d0)*&
1056 ! Bp_lm*R5p_lm)+((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*&
1057 ! (DPSUMo2*PI/6.d0)*R1p_lm) )*( (TH(M)*TH(LL))**(3.) )*&
1058 ! (NU_PM*GNUWDOTN(LL)-NU_PL*GNUWDOTN(M))
1059 !
1060 
1061  k_s_sum = k_s_sum + k_sm_lm
1062 ! Kth_s_sum = Kth_s_sum + Kth_sL_LM
1063 ! Kvel_s_sum = Kvel_s_sum + Kvel_s_LM
1064 ! Knu_s_sum = Knu_s_sum + Knu_s_LM
1065 
1066  ENDIF ! end if( LL .eq. M)/else
1067  ENDDO ! end do LL = 1, SMAX
1068 
1069 ! mth solids phase granular conductivity DOES NOT include particle mass
1070  k_1 = k_s_sum
1071 
1072 ! Note that the velocity term is not included here because it should
1073 ! become zero when dotted with the outward normal (i.e. no solids
1074 ! flux through the wall; assuming that the solids velocity at the
1075 ! wall in the normal direction is zero)
1076 ! CW = CW + Kvel_s_sum
1077 
1078 ! Note that the gradient in number density at the wall must be
1079 ! approximated with interior points since there is no value of
1080 ! number density associated with the wall location; the ghost
1081 ! cell values are undefined for volume fraction.
1082 ! CW = CW + Knu_s_sum
1083 
1084 ! The gradient in temperature of phase LL at the wall
1085 ! CW = CW + Kth_s_sum
1086 
1087 
1088  CASE (gd_1999)
1089 ! Note: k_boltz = M_PM
1090  d_pm = dp_avg(m)
1091  m_pm = (pi/6.d0)*(d_pm**3.)*ros_avg(m)
1092  nu_pm = (eps(m)*ros_avg(m))/m_pm
1093 
1094 ! Find pressure in the Mth solids phase
1095  press_star = one + 2.d0*(1.d0+c_e)*eps(m)*g0(m)
1096 
1097 ! find conductivity
1098  eta0 = 5.0d0*m_pm*dsqrt(th(m)/pi) / (16.d0*d_pm*d_pm)
1099 
1100  c_star = 32.0d0*(1.0d0 - c_e)*(1.d0 - 2.0d0*c_e*c_e) &
1101  / (81.d0 - 17.d0*c_e + 30.d0*c_e*c_e*(1.0d0-c_e))
1102 
1103  zeta0_star = (5.d0/12.d0)*g0(m)*(1.d0 - c_e*c_e) &
1104  * (1.d0 + (3.d0/32.d0)*c_star)
1105 
1106  kappa0 = (15.d0/4.d0)*eta0
1107 
1108  nu_kappa_star = (g0(m)/3.d0)*(1.d0+c_e) * ( 1.d0 + &
1109  (33.d0/16.d0)*(1.d0-c_e) + ((19.d0-3.d0*c_e)/1024.d0)*c_star)
1110 ! nu_mu_star = nu_kappa_star
1111 
1112  kappa_k_star = (2.d0/3.d0)*(1.d0 + 0.5d0*(1.d0+press_star)*c_star + &
1113  (3.d0/5.d0)*eps(m)*g0(m)*(1.d0+c_e)*(1.d0+c_e) * &
1114  (2.d0*c_e - 1.d0 + ( 0.5d0*(1.d0+c_e) - 5.d0/(3*(1.d0+c_e))) * &
1115  c_star ) ) / (nu_kappa_star - 2.d0*zeta0_star)
1116 
1117  kappa_star = kappa_k_star * (1.d0 + (6.d0/5.d0)*eps(m)* &
1118  g0(m)*(1.d0+c_e) ) + (256.d0/25.d0)*(eps(m)* &
1119  eps(m)/pi)*g0(m)*(1.d0+c_e)*(1.d0+(7.d0/32.d0)* &
1120  c_star)
1121 
1122  IF(switch == zero .OR. ro_g_avg == zero)THEN
1123  kgran_star = kappa0
1124  ELSEIF(th(m) .LT. small_number)THEN
1125  kgran_star = zero
1126  ELSE
1127  kgran_star = ros_avg(m)*eps(m)* g0(m)*th(m)* kappa0/ &
1128  (ros_avg(m)*g0(m)*eps(m)*th(m) + &
1129  1.2d0*dga*kappa0/ros_avg(m)) ! Note dgA is ~F_gs/ep_s
1130  ENDIF
1131 
1132 ! kgran_star includes particle mass
1133 ! mth solids phase granular conductivity includes particle mass
1134  k_1 = kgran_star * kappa_star
1135 
1136 ! transport coefficient of the Mth solids phase associated
1137 ! with gradient in volume fraction in heat flux
1138 ! qmu_k_star = 2.d0*( (1.d0+EPS(M)*DG_0DNU(EPS(M)))* &
1139 ! zeta0_star*kappa_k_star + ( (press_star/3.d0) + (2.d0/3.d0)* &
1140 ! EPS(M)*(1.d0+C_E) * (G0(M)+EPS(M)* &
1141 ! DG_0DNU(EPS(M))) )*c_star - (4.d0/5.d0)*EPS(M)* &
1142 ! G0(M)* (1.d0+(EPS(M)/2.d0)*DG_0DNU(EPS(M)))* &
1143 ! (1.d0+C_E) * ( C_E*(1.d0-C_E)+0.25d0*((4.d0/3.d0)+C_E* &
1144 ! (1.d0-C_E))*c_star ) ) / (2.d0*nu_kappa_star-3.d0*zeta0_star)
1145 ! qmu_star = qmu_k_star*(1.d0+(6.d0/5.d0)*EPS(M)*G0(M)*&
1146 ! (1.d0+C_E) )
1147 ! Kphis = (TH(M)*Kgran_star/NU_PM)*qmu_star
1148 
1149 ! defining granular pressure (for Jenkins BC)
1150  psotheta = ros_avg(m)*eps(m)*(one+2.d0*(one+c_e)*g0eps_avg)
1151  ! ~ROs_avg(m)*EPS(M)*press_star
1152  ps = ros_avg(m)*eps(m)*th(m)*(one+2.d0*(one+c_e)*g0eps_avg)
1153  ! ~ROs_avg(m)*EPS(M)*TH(M)*press_star
1154 
1155  CASE (gtsh_2012)
1156  d_pm = dp_avg(m)
1157  m_pm = (pi/6.d0)*(d_pm**3.)*ros_avg(m)
1158  nu_pm = (eps(m)*ros_avg(m))/m_pm
1159  chi = g0(m)
1160  vfrac = eps(m)
1161 
1162  re_g = epg*ro_g_avg*d_pm*vrel/mu_g_avg
1163  re_t = ro_g_avg*d_pm*dsqrt(th(m)) / mu_g_avg
1164  mu2_0 = dsqrt(2d0*pi) * chi * (one-c_e**2) ! eq. (6.22) GTSH theory
1165  mu4_0 = (4.5d0+c_e**2) * mu2_0 ! eq. (6.23) GTSH theory
1166  mu4_1 = (6.46875d0+0.9375d0*c_e**2)*mu2_0 + 2d0*dsqrt(2d0*pi)* &
1167  chi*(one+c_e) ! this is done to avoid /0 in case c_e = 1.0
1168  a2_gtshw = zero ! for eps = zero
1169  if((vfrac> small_number) .and. (th(m) > small_number)) then
1170  zeta_star = 4.5d0*dsqrt(2d0*pi)*(ro_g_avg/ros_avg(m))**2*re_g**2 * &
1171  s_star(vfrac,chi) / (vfrac*(one-vfrac)**2 * re_t**4)
1172  a2_gtshw = (5d0*mu2_0 - mu4_0) / (mu4_1 - 5d0* &
1173  (19d0/16d0*mu2_0 - 1.5d0*zeta_star))
1174  endif
1175  eta0 = 0.3125d0/(dsqrt(pi)*d_pm**2)*m_pm*dsqrt(th(m))
1176  nu0 = (96.d0/5.d0)*(vfrac/d_pm)*dsqrt(th(m)/pi)
1177  nuk = nu0*(one+c_e)/3d0*chi*( one+2.0625d0*(one-c_e)+ &
1178  ((947d0-579*c_e)/256d0*a2_gtshw) )
1179  kth0 = 3.75d0*eta0/m_pm
1180  edt_s = 4d0/3d0*dsqrt(pi)*(one-c_e**2)*chi* &
1181  (one+0.1875d0*a2_gtshw)*nu_pm*d_pm**2*dsqrt(th(m))
1182  kthk = zero
1183  if(vfrac > small_number) kthk = 2d0/3d0*kth0*nu0/(nuk - 2d0*edt_s) * &
1184  (one+2d0*a2_gtshw+0.6d0*vfrac*chi* &
1185  (one+c_e)**2*(2*c_e-one+a2_gtshw*(one+c_e)))
1186 
1187 ! defining conductivity K_1 at walls equivalent to Kth_s(IJK,M) in calc_mu_s
1188 ! mth solids phase granular conductivity DOES NOT include particle mass
1189  k_1 = kthk*(one+1.2d0*vfrac*chi*(one+c_e)) + (10.24d0/pi* &
1190  vfrac**2*chi*(one+c_e)*(one+0.4375d0*a2_gtshw)*kth0)
1191 
1192 ! defining granular pressure (for Jenkins BC)
1193  psotheta = ros_avg(m)*eps(m)*(one+2.d0*(one+c_e)*g0eps_avg)
1194  ps = ros_avg(m)*eps(m)*th(m)*(one+2.d0*(one+c_e)*g0eps_avg)
1195 
1196 
1197  CASE DEFAULT
1198 ! should never hit this
1199  WRITE (*, '(A)') 'BC_THETA => THETA_HW_CW'
1200  WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', kt_type
1201  call mfix_exit(mype)
1202  END SELECT
1203 
1204 
1205 ! setting the coefficients for JJ BC
1206 ! -------------------------------------------------------------------
1207 ! GW = granular conductivity
1208 ! Hw = dissipation due to particle-wall collision
1209 ! Cw = generation due to particle-wall slip
1210  SELECT CASE (kt_type_enum)
1211  CASE (lun_1984, simonin_1996, ahmadi_1995, gd_1999)
1212 ! KTs without particle mass in their definition of granular temperature
1213 ! and with mass in their conductivity
1214 ! and theta = granular temperature
1215  gw = k_1
1216 
1217 ! Jenkins BC implemented for these KT models
1218  IF(jenkins) THEN
1219  hw = 3.d0/8.d0*dsqrt(3.d0*th(m))*((1d0-e_w))*&
1220  psotheta
1221 
1222 ! As I understand from soil mechanic papers, the coefficient mu in
1223 ! Jenkins paper is tan_Phi_w. T. See for example, G.I. Tardos, PT,
1224 ! 92 (1997), 61-74, equation (1). sof
1225  cw = tan_phi_w*tan_phi_w*(one+e_w)*21.d0/16.d0*&
1226  dsqrt(3.d0*th(m)) * ps
1227 
1228  ELSE
1229  hw = (pi*dsqrt(3d0)/(4.d0*(one-ep_star_avg)))*&
1230  (1d0-e_w*e_w)*&
1231  ros_avg(m)*eps(m)*g0(m)*dsqrt(th(m))
1232 
1233  IF(.NOT. bc_jj_m) THEN
1234  cw = (pi*dsqrt(3d0)/(6.d0*(one-ep_star_avg)))*phip*&
1235  ros_avg(m)*eps(m)*g0(m)*dsqrt(th(m))*vslipsq
1236  ELSE
1237  cw = (pi*dsqrt(3d0)/(6.d0*(one-ep_star_avg)))*&
1238  phip_jj(dsqrt(vslipsq),th(m))*ros_avg(m)*&
1239  eps(m)*g0(m)*dsqrt(th(m))*vslipsq
1240  ENDIF
1241  ENDIF ! end if(Jenkins)/else
1242 
1243 
1244  CASE (gtsh_2012)
1245 ! KTs with particle mass in their definition of granular temperature
1246 ! and without mass in their conductivity
1247 ! and theta = granular temperature/mass
1248  gw = m_pm * k_1
1249 
1250  hw = (pi*dsqrt(3d0)/(4.d0*(one-ep_star_avg)))*&
1251  (1d0-e_w*e_w)*&
1252  ros_avg(m)*eps(m)*g0(m)*dsqrt(th(m))
1253 
1254  IF(.NOT. bc_jj_m) THEN
1255  cw = (pi*dsqrt(3d0)/(6.d0*(one-ep_star_avg)))*phip*&
1256  ros_avg(m)*eps(m)*g0(m)*dsqrt(th(m))*vslipsq
1257  ELSE
1258  cw = (pi*dsqrt(3d0)/(6.d0*(one-ep_star_avg)))*&
1259  phip_jj(dsqrt(vslipsq),th(m))*ros_avg(m)*&
1260  eps(m)*g0(m)*dsqrt(th(m))*vslipsq
1261  ENDIF
1262 
1263 
1264  CASE (ia_2005)
1265 ! KTs with particle mass in their definition of granular temperature
1266 ! and without mass in their conductivity
1267 ! and theta = granular temperature
1268  gw = m_pm * k_1
1269 
1270  hw = (pi*dsqrt(3.d0)/(4.d0*(one-ep_star_avg)))*&
1271  (1.d0-e_w*e_w)*&
1272  ros_avg(m)*eps(m)*g0(m)*dsqrt(th(m)/m_pm)
1273 
1274  IF(.NOT. bc_jj_m) THEN
1275  cw = (pi*dsqrt(3.d0)/(6.d0*(one-ep_star_avg)))*phip*&
1276  ros_avg(m)*eps(m)*g0(m)*dsqrt(th(m)*m_pm)*vslipsq
1277  ELSE
1278  cw = (pi*dsqrt(3.d0)/(6.d0*(one-ep_star_avg)))*&
1279  phip_jj(dsqrt(vslipsq),th(m))*ros_avg(m)*&
1280  eps(m)*g0(m)*dsqrt(th(m)*m_pm)*vslipsq
1281  ENDIF
1282 
1283  CASE DEFAULT
1284 ! should never hit this
1285  WRITE (*, '(A)') 'BC_THETA => THETA_HW_CW'
1286  WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', kt_type
1287  call mfix_exit(mype)
1288  END SELECT
1289 
1290  RETURN
1291  END SUBROUTINE theta_hw_cw
1292 
1293 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1294 ! C
1295 ! Function: PHIP_JJ C
1296 ! Purpose: Calculate the specularity coefficient for JJ BC C
1297 ! C
1298 ! Author: Tingwen Li C
1299 ! C
1300 ! C
1301 ! Literature/Document References: C
1302 ! Li, T., and Benyahia, S., Revisiting Johnson and Jackson C
1303 ! boundary conditions for granular flows, AIChE Journal, 2012, C
1304 ! Vol 58 (7), pp. 2058-2068 C
1305 ! C
1306 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1307 
1308 
1309  DOUBLE PRECISION FUNCTION phip_jj(uslip,g_theta)
1310 !-----------------------------------------------
1311 ! Modules
1312 !-----------------------------------------------
1313  use constant, only : pi, k4phi, phip0
1314  implicit none
1315 !-----------------------------------------------
1316 ! Dummy arguments
1317 !-----------------------------------------------
1318  double precision, INTENT(IN) :: uslip
1319  double precision, INTENT(IN) :: g_theta
1320 !-----------------------------------------------
1321 ! Local variables
1322 !-----------------------------------------------
1323  double precision :: r4phi
1324 !-----------------------------------------------
1325 
1326  !k4phi=7.d0/2.d0*mu4phi*(1.d0+e_w)
1327  r4phi=uslip/dsqrt(3.d0*g_theta)
1328 
1329  IF(r4phi .le. 4.d0*k4phi/7.d0/dsqrt(6.d0*pi)/phip0) THEN
1330  phip_jj=-7.d0*dsqrt(6.d0*pi)*phip0**2/8.d0/k4phi*r4phi+phip0
1331  ELSE
1332  phip_jj=2.d0/7.d0*k4phi/r4phi/dsqrt(6.d0*pi)
1333  ENDIF
1334 
1335  RETURN
1336 
1337  END FUNCTION phip_jj
bc::bc_ww_s
double precision, dimension(dimension_bc, dim_m) bc_ww_s
Definition: bc_mod.f:328
bc::bc_k_b
integer, dimension(dimension_bc) bc_k_b
Definition: bc_mod.f:70
mpi_utility
Definition: mpi_utility_mod.f:6
constant::c_e
double precision c_e
Definition: constant_mod.f:105
constant::e_w
double precision e_w
Definition: constant_mod.f:85
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
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
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
write_error
subroutine write_error(NAME, LINE, LMAX)
Definition: write_error.f:21
bc::bc_i_w
integer, dimension(dimension_bc) bc_i_w
Definition: bc_mod.f:54
fldvar::w_s
double precision, dimension(:,:), allocatable w_s
Definition: fldvar_mod.f:117
bc::bc_j_n
integer, dimension(dimension_bc) bc_j_n
Definition: bc_mod.f:66
rxns
Definition: rxns_mod.f:1
indices::im1
integer, dimension(:), allocatable im1
Definition: indices_mod.f:50
constant::switch
double precision, parameter switch
Definition: constant_mod.f:53
param::dimension_bc
integer, parameter dimension_bc
Definition: param_mod.f:61
bc::bc_type_enum
integer, dimension(dimension_bc) bc_type_enum
Definition: bc_mod.f:146
turb::k_12
double precision, dimension(:), allocatable k_12
Definition: turb_mod.f:17
rdf::g_0
double precision function g_0(IJK, M1, M2)
Definition: rdf_mod.f:240
param::east
integer east
Definition: param_mod.f:29
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
run::phip_out_iter
integer phip_out_iter
Definition: run_mod.f:172
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
compar::pe_io
integer pe_io
Definition: compar_mod.f:30
geometry::ody_n
double precision, dimension(:), allocatable ody_n
Definition: geometry_mod.f:123
fun_avg
Definition: fun_avg_mod.f:1
run::jenkins
logical jenkins
Definition: run_mod.f:166
theta_hw_cw
subroutine theta_hw_cw(g0, EPS, EPG, ep_star_avg, g0EPs_avg, TH, Mu_g_avg, RO_g_avg, ROs_avg, DP
Definition: bc_theta.f:743
fldvar::d_p
double precision, dimension(:,:), allocatable d_p
Definition: fldvar_mod.f:57
constant::k4phi
double precision k4phi
Definition: constant_mod.f:83
exit::mfix_exit
subroutine mfix_exit(myID, normal_termination)
Definition: exit.f:5
bc::bc_k_t
integer, dimension(dimension_bc) bc_k_t
Definition: bc_mod.f:74
compar::dead_cell_at
logical, dimension(:,:,:), allocatable dead_cell_at
Definition: compar_mod.f:127
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
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
bc::bc_j_s
integer, dimension(dimension_bc) bc_j_s
Definition: bc_mod.f:62
exit
Definition: exit.f:2
geometry::ox
double precision, dimension(:), allocatable ox
Definition: geometry_mod.f:140
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
param::north
integer north
Definition: param_mod.f:37
fldvar::v_g
double precision, dimension(:), allocatable v_g
Definition: fldvar_mod.f:99
constant::phip
double precision phip
Definition: constant_mod.f:79
bc::bc_defined
logical, dimension(dimension_bc) bc_defined
Definition: bc_mod.f:207
calc_theta_bc
subroutine calc_theta_bc(IJK1, IJK2, FCELL, M, L, Gw, Hw, Cw)
Definition: bc_theta.f:292
param::south
integer south
Definition: param_mod.f:41
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
toleranc::dil_ep_s
double precision, parameter dil_ep_s
Definition: toleranc_mod.f:24
fldvar::ro_s
double precision, dimension(:,:), allocatable ro_s
Definition: fldvar_mod.f:45
indices::km1
integer, dimension(:), allocatable km1
Definition: indices_mod.f:52
compar::mype
integer mype
Definition: compar_mod.f:24
rxns::reactionrates
double precision, dimension(:,:), allocatable reactionrates
Definition: rxns_mod.f:7
physprop::mu_g
double precision, dimension(:), allocatable mu_g
Definition: physprop_mod.f:68
param::west
integer west
Definition: param_mod.f:33
run::phip_out_jj
logical phip_out_jj
Definition: run_mod.f:170
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
param::top
integer top
Definition: param_mod.f:45
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
bc::bc_jj_ps
integer, dimension(dimension_bc) bc_jj_ps
Definition: bc_mod.f:212
phip_jj
double precision function phip_jj(uslip, g_theta)
Definition: bc_theta.f:1310
mpi_utility::exitmpi
subroutine exitmpi(myid)
Definition: mpi_utility_mod.f:5437
constant::pi
double precision, parameter pi
Definition: constant_mod.f:158
constant::phip0
double precision phip0
Definition: constant_mod.f:81
param::dimension_m
integer dimension_m
Definition: param_mod.f:18
geometry::odz_t
double precision, dimension(:), allocatable odz_t
Definition: geometry_mod.f:125
bc::bc_vw_s
double precision, dimension(dimension_bc, dim_m) bc_vw_s
Definition: bc_mod.f:325
fldvar::ro_g
double precision, dimension(:), allocatable ro_g
Definition: fldvar_mod.f:32
param::bottom
integer bottom
Definition: param_mod.f:49
bc::bc_i_e
integer, dimension(dimension_bc) bc_i_e
Definition: bc_mod.f:58
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
bc_theta
subroutine bc_theta(M, A_m, B_m)
Definition: bc_theta.f:31
bc
Definition: bc_mod.f:23
constant::switch_ia
logical, parameter switch_ia
Definition: constant_mod.f:74