File: RELATIVE:/../../../mfix.git/model/bc_theta.f

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 matrix
39           USE scales
40           USE constant
41           USE toleranc
42           USE run
43           USE physprop
44           USE fldvar
45           USE visc_s
46           USE geometry
47           USE output
48           USE indices
49           USE bc
50           USE compar
51           USE mpi_utility
52           USE fun_avg
53           USE functions
54           IMPLICIT NONE
55     !-----------------------------------------------
56     ! Dummy arguments
57     !-----------------------------------------------
58     ! Solids phase index
59           INTEGER :: M
60     ! Septadiagonal matrix A_m
61           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
62     ! Vector b_m
63           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
64     !-----------------------------------------------
65     ! Local variables
66     !-----------------------------------------------
67     ! Boundary condition
68           INTEGER :: L
69     ! Indices
70           INTEGER :: I,  J, K, I1, I2, J1, J2, K1, K2, IJK,&
71                      IM, JM, KM
72     ! Granular energy coefficient
73           DOUBLE PRECISION :: Gw, Hw, Cw
74     !-----------------------------------------------
75     
76     ! Setup Johnson and Jackson Pseudo-thermal temp B.C.
77           DO L = 1, DIMENSION_BC
78             IF( BC_DEFINED(L) ) THEN
79               IF(BC_TYPE(L) .EQ. 'NO_SLIP_WALL' .OR.&
80                  BC_TYPE(L) .EQ. 'FREE_SLIP_WALL' .OR.&
81                  BC_TYPE(L) .EQ. 'PAR_SLIP_WALL' ) THEN
82                 I1 = BC_I_w(L)
83                 I2 = BC_I_e(L)
84                 J1 = BC_J_s(L)
85                 J2 = BC_J_n(L)
86                 K1 = BC_K_b(L)
87                 K2 = BC_K_t(L)
88     
89                 IF(BC_JJ_PS(L).GT.0) THEN
90                   DO K = K1, K2
91                   DO J = J1, J2
92                   DO I = I1, I2
93     
94                     IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
95                     IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
96                     IJK   = FUNIJK(I, J, K)
97                     IF (FLOW_AT(IJK)) CYCLE !checks for pressure outlets
98                     IM    = Im1(I)
99                     JM    = Jm1(J)
100                     KM    = Km1(K)
101     ! Setting the temperature to zero - recall the center coefficient
102     ! and source vector are negative. The off-diagonal coefficients
103     ! are positive.
104                     A_m(IJK, e, M) =  ZERO
105                     A_m(IJK, w, M) =  ZERO
106                     A_m(IJK, n, M) =  ZERO
107                     A_m(IJK, s, M) =  ZERO
108                     A_m(IJK, t, M) =  ZERO
109                     A_m(IJK, b, M) =  ZERO
110                     A_m(IJK, 0, M) = -ONE
111                     b_m(IJK, M)    =  ZERO
112     
113     ! Checking if the west wall
114                     IF(FLUID_AT(EAST_OF(IJK)))THEN
115                        IF(EP_s(EAST_OF(IJK),M).LE.DIL_EP_s)THEN
116                           A_m(IJK, e, M) = ONE
117                        ELSE
118                           IF (BC_JJ_PS(L).EQ.3) THEN
119     ! Setting the wall temperature equal to the adjacent fluid cell
120     ! temperature (i.e. zero flux)
121                              Gw = 1d0
122                              Hw = 0d0
123                              Cw = 0d0
124                           ELSE
125     ! Setting the wall temperature based on Johnson and Jackson
126     ! Pseudo-thermal temp B.C.
127                              CALL CALC_THETA_BC(IJK,EAST_OF(IJK),'E',&
128                                 M,L,gw,hw,cw)
129     ! Overwriting calculated value of cw
130                             IF (BC_JJ_PS(L).EQ.2) cw=0d0
131                           ENDIF
132     
133                           A_m(IJK, e, M) = -(HALF*hw - oDX_E(I)*gw)
134                           A_m(IJK, 0, M) = -(HALF*hw + oDX_E(I)*gw)
135     
136                           IF (BC_JJ_PS(L) .EQ. 1) THEN
137                             b_m(IJK, M) = -cw
138                           ELSE
139                             b_m(IJK,M) = ZERO
140                           ENDIF
141                        ENDIF
142     
143                     ELSEIF(FLUID_AT(WEST_OF(IJK)))THEN
144                        IF(EP_s(WEST_OF(IJK),M).LE.DIL_EP_s)THEN
145                           A_m(IJK, w, M) = ONE
146                        ELSE
147                          IF (BC_JJ_PS(L).EQ.3) THEN
148                             Gw = 1d0
149                             Hw = 0d0
150                             Cw = 0d0
151                          ELSE
152                             CALL CALC_THETA_BC(IJK,WEST_OF(IJK),'W',&
153                                M,L,gw,hw,cw)
154                             IF (BC_JJ_PS(L).EQ.2) cw=0d0
155                          ENDIF
156     
157                          A_m(IJK, w, M) = -(HALF*hw - oDX_E(IM)*gw)
158                          A_m(IJK, 0, M) = -(HALF*hw + oDX_E(IM)*gw)
159     
160                          IF (BC_JJ_PS(L) .EQ. 1) THEN
161                             b_m(IJK, M) = -cw
162                          ELSE
163                             b_m(IJK,M) = ZERO
164                          ENDIF
165                        ENDIF
166     
167                     ELSEIF(FLUID_AT(NORTH_OF(IJK)))THEN
168                        IF(EP_s(NORTH_OF(IJK),M).LE.DIL_EP_s)THEN
169                           A_m(IJK, n, M) = ONE
170                        ELSE
171                          IF (BC_JJ_PS(L).EQ.3) THEN
172                             Gw = 1d0
173                             Hw = 0d0
174                             Cw = 0d0
175                          ELSE
176                             CALL CALC_THETA_BC(IJK,NORTH_OF(IJK),'N',&
177                                M,L,gw,hw,cw)
178                             IF (BC_JJ_PS(L).EQ.2) cw=0d0
179                          ENDIF
180     
181                          A_m(IJK, n, M) = -(HALF*hw - oDY_N(J)*gw)
182                          A_m(IJK, 0, M) = -(HALF*hw + oDY_N(J)*gw)
183     
184                          IF (BC_JJ_PS(L) .EQ. 1) THEN
185                             b_m(IJK, M) = -cw
186                          ELSE
187                             b_m(IJK,M) = ZERO
188                          ENDIF
189                        ENDIF
190     
191                     ELSEIF(FLUID_AT(SOUTH_OF(IJK)))THEN
192                        IF(EP_s(SOUTH_OF(IJK),M).LE.DIL_EP_s)THEN
193                           A_m(IJK, s, M) = ONE
194                        ELSE
195                          IF (BC_JJ_PS(L).EQ.3) THEN
196                             Gw = 1d0
197                             Hw = 0d0
198                             Cw = 0d0
199                          ELSE
200                             CALL CALC_THETA_BC(IJK,SOUTH_OF(IJK),'S',&
201                                M,L,gw,hw,cw)
202                             IF (BC_JJ_PS(L).EQ.2) cw=0d0
203                          ENDIF
204     
205                          A_m(IJK, s, M) = -(HALF*hw - oDY_N(JM)*gw)
206                          A_m(IJK, 0, M) = -(HALF*hw + oDY_N(JM)*gw)
207     
208                          IF (BC_JJ_PS(L) .EQ. 1) THEN
209                             b_m(IJK, M) = -cw
210                          ELSE
211                             b_m(IJK,M) = ZERO
212                          ENDIF
213                        ENDIF
214     
215                     ELSEIF(FLUID_AT(TOP_OF(IJK)))THEN
216                        IF(EP_s(TOP_OF(IJK),M).LE.DIL_EP_s)THEN
217                           A_m(IJK, t, M) = ONE
218                        ELSE
219                          IF (BC_JJ_PS(L).EQ.3) THEN
220                             Gw = 1d0
221                             Hw = 0d0
222                             Cw = 0d0
223                          ELSE
224                             CALL CALC_THETA_BC(IJK,TOP_OF(IJK),'T',&
225                                M,L,gw,hw,cw)
226                             IF (BC_JJ_PS(L).EQ.2) cw=0d0
227                          ENDIF
228     
229                          A_m(IJK, t, M) = -(HALF*hw - oX(I)*oDZ_T(K)*gw)
230                          A_m(IJK, 0, M) = -(HALF*hw + oX(I)*oDZ_T(K)*gw)
231     
232                          IF (BC_JJ_PS(L) .EQ. 1) THEN
233                             b_m(IJK, M) = -cw
234                          ELSE
235                             b_m(IJK,M) = ZERO
236                          ENDIF
237                        ENDIF
238     
239                     ELSEIF(FLUID_AT(BOTTOM_OF(IJK)))THEN
240                        IF(EP_s(BOTTOM_OF(IJK),M).LE.DIL_EP_s)THEN
241                           A_m(IJK, b , M) = ONE
242                        ELSE
243                          IF (BC_JJ_PS(L).EQ.3) THEN
244                             Gw = 1d0
245                             Hw = 0d0
246                             Cw = 0d0
247                          ELSE
248                             CALL CALC_THETA_BC(IJK,BOTTOM_OF(IJK),'B',&
249                                M,L,gw,hw,cw)
250                             IF (BC_JJ_PS(L).EQ.2) cw=0d0
251                          ENDIF
252     
253                          A_m(IJK, b, M) = -(HALF*hw - oX(I)*oDZ_T(KM)*gw)
254                          A_m(IJK, 0, M) = -(HALF*hw + oX(I)*oDZ_T(KM)*gw)
255     
256                          IF (BC_JJ_PS(L) .EQ. 1) THEN
257                             b_m(IJK, M) = -cw
258                          ELSE
259                             b_m(IJK,M) = ZERO
260                          ENDIF
261                        ENDIF
262                     ENDIF
263     
264                   ENDDO   ! end do over i
265                   ENDDO   ! end do over j
266                   ENDDO   ! end do over k
267     
268                 ENDIF   ! end if (bc_jj_ps>0)
269               ENDIF  ! end if nsw, fsw or psw
270             ENDIF   ! end if (bc_defined)
271           ENDDO   ! end L do loop over dimension_bc
272     
273           RETURN
274           END SUBROUTINE BC_THETA
275     
276     
277     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
278     !                                                                      C
279     !  Subroutine: CALC_THETA_BC                                           C
280     !                                                                      C
281     !  Purpose: Implementation of Johnson & Jackson boundary conditions    C
282     !  for the pseudo-thermal temperature.                                 C
283     !                                                                      C
284     !                                                                      C
285     !  Author: Kapil Agrawal, Princeton University        Date: 14-MAR-98  C
286     !  Reviewer:                                          Date:            C
287     !                                                                      C
288     !                                                                      C
289     !                                                                      C
290     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
291     
292           SUBROUTINE CALC_THETA_BC(IJK1,IJK2,FCELL,M,L,Gw,Hw,Cw)
293     
294     !-----------------------------------------------
295     ! Modules
296     !-----------------------------------------------
297           USE bc
298           USE compar
299           USE constant
300           USE fldvar
301           USE fun_avg
302           USE functions
303           USE geometry
304           USE indices
305           USE mpi_utility
306           USE param
307           USE param1
308           USE physprop
309           USE rdf
310           USE run
311           USE rxns
312           USE toleranc
313           USE turb
314           USE visc_s
315           IMPLICIT NONE
316     !-----------------------------------------------
317     ! Dummy Arguments
318     !-----------------------------------------------
319     ! IJK indices for wall cell (ijk1) and fluid cell (ijk2)
320           INTEGER, INTENT(IN) :: IJK1, IJK2
321     ! The location (e,w,n...) of fluid cell
322           CHARACTER, INTENT(IN) :: FCELL
323     ! Solids phase index
324           INTEGER, INTENT(IN) :: M
325     ! Index corresponding to boundary condition
326           INTEGER, INTENT(IN) ::  L
327     ! Granular energy coefficient
328           DOUBLE PRECISION, INTENT(INOUT) :: Gw, Hw, Cw
329     !-----------------------------------------------
330     ! Local Variables
331     !-----------------------------------------------
332     ! IJK indices for cells
333           INTEGER ::  IJK3
334     ! Solids phase index
335           INTEGER :: MM
336     ! Average scalars
337           DOUBLE PRECISION :: EPg_avg, Mu_g_avg, RO_g_avg, smallTheta
338     ! void fraction at packing
339           DOUBLE PRECISION :: ep_star_avg
340     ! Average scalars modified to include all solid phases
341           DOUBLE PRECISION :: EPs_avg(DIMENSION_M), DP_avg(DIMENSION_M), &
342                               TH_avg(DIMENSION_M)
343     ! Average density of solids phase M
344           DOUBLE PRECISION ROs_avg(DIMENSION_M)
345     ! Average Simonin variables
346           DOUBLE PRECISION :: K_12_avg, Tau_12_avg
347     ! values of U_sm, V_sm, W_sm at appropriate place on boundary wall
348           DOUBLE PRECISION :: USCM, VSCM, WSCM
349     ! values of U_g, V_g, W_g at appropriate place on boundary wall
350           DOUBLE PRECISION :: UGC, VGC, WGC
351     ! Magnitude of gas-solids relative velocity
352           DOUBLE PRECISION :: VREL
353     ! Square of slip velocity between wall and particles
354           DOUBLE PRECISION :: VSLIPSQ
355     ! Wall velocity dot outward normal for all solids phases
356           DOUBLE PRECISION :: VWDOTN (DIMENSION_M)
357     ! Gradient in number density at the wall dot with outward normal for
358     ! all solids phases
359           DOUBLE PRECISION :: GNUWDOTN (DIMENSION_M)
360     ! Gradient in granular temperature at the wall dot with outward normal
361     ! for all solids phases
362           DOUBLE PRECISION :: GTWDOTN (DIMENSION_M)
363     ! Error message
364           CHARACTER(LEN=80)     LINE(1)
365     ! Radial distribution function
366           DOUBLE PRECISION :: g0(DIMENSION_M)
367     ! Sum of eps*G_0
368           DOUBLE PRECISION :: g0EPs_avg
369     !-----------------------------------------------
370     ! External functions
371     !-----------------------------------------------
372     ! Variable specularity coefficient
373           DOUBLE PRECISION :: PHIP_JJ
374     !-----------------------------------------------
375     
376     ! Note: EP_s, MU_g, and RO_g are undefined at IJK1 (wall cell).
377     !       Hence IJK2 (fluid cell) is used in averages.
378     
379           smallTheta = (to_SI)**4 * ZERO_EP_S
380     
381     ! set location independent quantities
382           EPg_avg = EP_g(IJK2)
383           ep_star_avg = EP_star_array(IJK2)
384           Mu_g_avg = Mu_g(IJK2)
385           RO_g_avg = RO_g(IJK2)
386           g0EPs_avg = ZERO
387     
388     ! added for Simonin model (sof)
389           IF(KT_TYPE_ENUM == SIMONIN_1996) THEN
390              K_12_avg = K_12(IJK2)
391              Tau_12_avg = Tau_12(IJK2)
392           ELSE
393              K_12_avg = ZERO
394              Tau_12_avg = ZERO
395           ENDIF
396     
397           DO MM = 1, SMAX
398              g0(MM)      = G_0(IJK2, M, MM)
399              EPs_avg(MM) = EP_s(IJK2,MM)
400              DP_avg(MM)  = D_P(IJK2,MM)
401              g0EPs_avg   = g0EPs_avg + G_0(IJK2, M, MM)*EP_s(IJK2,MM)
402              ROs_avg(MM) = RO_S(IJK2,MM)
403           ENDDO
404     
405     
406           IF(FCELL .EQ. 'N')THEN
407              DO MM = 1, SMAX
408                 TH_avg(MM) = AVG_Y(Theta_m(IJK1,MM),Theta_m(IJK2,MM),J_OF(IJK1))
409                 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta ! for some corner cells
410     
411     ! added for IA (2005) theory:
412     ! include -1 since normal vector is pointing south (-y)
413     ! velocity at wall (i,j+1/2,k relative to ijk1) dot with the normal
414     ! vector at the wall
415                 VWDOTN(MM) = -1.d0*V_S(IJK1,MM)
416     
417     ! gradient in number density at wall (i,j+1/2,k relative to ijk1) dot
418     ! with the normal vector at the wall. since nu is undefined at ijk1,
419     ! approximate gradient with interior points: ijk2 and i,j+1,k relative
420     ! to ijk2
421                 IJK3 = NORTH_OF(IJK2)
422                 GNUWDOTN(MM) = -1.d0*(6.d0/(PI*DP_avg(MM)))*&
423                      (EP_s(IJK3,MM)-EP_s(IJK2,MM))*oDY_N(J_OF(IJK2))
424     
425     ! gradient in granular temperature at wall (i,j+1/2,k relative to ijk1) dot
426     ! with the normal vector at the wall.
427                 GTWDOTN(MM) = -1.d0*(Theta_m(IJK3,MM)-Theta_m(IJK2,MM))*&
428                      oDY_N(J_OF(IJK2))
429              ENDDO
430     
431     ! Calculate velocity components at i, j+1/2, k (relative to IJK1)
432              UGC  = AVG_Y(AVG_X_E(U_g(IM_OF(IJK1)),U_g(IJK1),I_OF(IJK1)),&
433                          AVG_X_E(U_g(IM_OF(IJK2)),U_g(IJK2),I_OF(IJK2)),&
434                          J_OF(IJK1))
435              VGC  = V_g(IJK1)
436              WGC  = AVG_Y(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
437                          AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
438                          J_OF(IJK1))
439              USCM = AVG_Y(AVG_X_E(U_s(IM_OF(IJK1),M),U_s(IJK1,M),I_OF(IJK1)),&
440                          AVG_X_E(U_s(IM_OF(IJK2),M),U_s(IJK2,M),I_OF(IJK2)),&
441                          J_OF(IJK1))
442              VSCM = V_s(IJK1,M)
443              WSCM = AVG_Y(AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
444                          AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
445                          J_OF(IJK1))
446     
447     ! slip velocity at the wall
448              VSLIPSQ=(WSCM-BC_Ww_s(L, M))**2 + (USCM-BC_Uw_s(L, M))**2
449     
450     
451           ELSEIF(FCELL .EQ. 'S')THEN
452              DO MM = 1, SMAX
453                 TH_avg(MM) = AVG_Y(Theta_m(IJK2, MM),Theta_m(IJK1, MM),J_OF(IJK2))
454                 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta ! for some corner cells
455     
456     ! added for IA (2005) theory:
457     ! include 1 since normal vector is pointing north (+y)
458     ! velocity at wall (i,j+1/2,k relative to ijk2) dot with the normal
459     ! vector at the wall.
460                 VWDOTN(MM) = 1.d0*V_S(IJK2,MM)
461     
462     ! gradient in number density at wall (i,j+1/2,k relative to ijk2) dot
463     ! with the normal vector at the wall. since nu is undefined at ijk1,
464     ! approximate gradient with interior points: ijk2 and i,j-1,k relative
465     ! to ijk2
466                 IJK3 = SOUTH_OF(IJK2)
467                 GNUWDOTN(MM) = 1.d0*(6.d0/(PI*DP_avg(MM)))*&
468                      (EP_s(IJK2,MM)-EP_s(IJK3,MM))*oDY_N(J_OF(IJK3))
469     
470     ! gradient in granular temperature at wall (i,j+1/2,k relative to ijk2)
471     ! dot with the normal vector at the wall.
472                 GTWDOTN(MM) = 1.d0*(Theta_m(IJK2,MM)-Theta_m(IJK3,MM))*&
473                      oDY_N(J_OF(IJK3))
474              ENDDO
475     
476     ! Calculate velocity components at i, j+1/2, k (relative to IJK2)
477              UGC  = AVG_Y(AVG_X_E(U_g(IM_OF(IJK2)),U_g(IJK2),I_OF(IJK2)),&
478                          AVG_X_E(U_g(IM_OF(IJK1)),U_g(IJK1),I_OF(IJK1)),&
479                          J_OF(IJK2))
480              VGC  = V_g(IJK2)
481              WGC  = AVG_Y(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
482                          AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
483                          J_OF(IJK2))
484              USCM = AVG_Y(AVG_X_E(U_s(IM_OF(IJK2),M),U_s(IJK2,M),I_OF(IJK2)),&
485                          AVG_X_E(U_s(IM_OF(IJK1),M),U_s(IJK1,M),I_OF(IJK1)),&
486                          J_OF(IJK2))
487              VSCM = V_s(IJK2,M)
488              WSCM = AVG_Y(AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
489                          AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
490                          J_OF(IJK2))
491     
492     ! slip velocity at the wall
493              VSLIPSQ=(WSCM-BC_Ww_s(L, M))**2 + (USCM-BC_Uw_s(L, M))**2
494     
495     
496           ELSEIF(FCELL== 'E')THEN
497               DO MM = 1, SMAX
498                 TH_avg(MM) = AVG_X(Theta_m(IJK1,MM),Theta_m(IJK2,MM),I_OF(IJK1))
499                 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta ! for some corner cells
500     
501     ! added for IA (2005) theory:
502     ! include -1 since normal vector is pointing west (-x)
503     ! velocity at wall (i+1/2,j,k relative to ijk1) dot with the normal
504     ! vector at the wall.
505                 VWDOTN(MM) = -1.d0*U_S(IJK1,MM)
506     
507     ! gradient in number density at wall (i+1/2,j,k relative to ijk1) dot
508     ! with the normal vector at the wall. since nu is undefined at ijk1,
509     ! approximate gradient with interior points: ijk2 and i,i+1,k relative
510     ! to ijk2
511                 IJK3 = EAST_OF(IJK2)
512                 GNUWDOTN(MM) = -1.d0*(6.d0/(PI*DP_avg(MM)))*&
513                      (EP_s(IJK3,MM)-EP_s(IJK2,MM))*oDX_E(I_OF(IJK2))
514     
515     ! gradient in granular temperature at wall (i+1/2,j,k relative to ijk1)
516     ! dot with the normal vector at the wall.
517                 GTWDOTN(MM) = -1.d0*(Theta_m(IJK3,MM)-Theta_m(IJK2,MM))*&
518                      oDX_E(I_OF(IJK2))
519              ENDDO
520     
521     ! Calculate velocity components at i+1/2, j, k (relative to IJK1)
522              UGC  = U_g(IJK1)
523              VGC  = AVG_X(AVG_Y_N(V_g(JM_OF(IJK1)), V_g(IJK1)),&
524                          AVG_Y_N(V_g(JM_OF(IJK2)), V_g(IJK2)),&
525                          I_OF(IJK1))
526              WGC  = AVG_X(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
527                          AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
528                          I_OF(IJK1))
529              USCM = U_s(IJK1,M)
530              VSCM = AVG_X(AVG_Y_N(V_s(JM_OF(IJK1),M), V_s(IJK1,M)),&
531                          AVG_Y_N(V_s(JM_OF(IJK2),M), V_s(IJK2,M)),&
532                          I_OF(IJK1))
533              WSCM = AVG_X(AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
534                          AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
535                          I_OF(IJK1))
536     
537     ! slip velocity at the wall
538              VSLIPSQ=(WSCM-BC_Ww_s(L, M))**2 + (VSCM-BC_Vw_s(L, M))**2
539     
540     
541           ELSEIF(FCELL== 'W')THEN
542               DO MM = 1, SMAX
543                 TH_avg(MM) = AVG_X(Theta_m(IJK2,MM),Theta_m(IJK1,MM),I_OF(IJK2))
544                 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta ! for some corner cells
545     
546     ! added for IA (2005) theory:
547     ! include 1 since normal vector is pointing west (+x)
548     ! velocity at wall (i+1/2,j,k relative to ijk2) dot with the normal
549     ! vector at the wall.
550                 VWDOTN(MM) = 1.d0*U_S(IJK2,MM)
551     
552     ! gradient in number density at wall (i+1/2,j,k relative to ijk2) dot
553     ! with the normal vector at the wall. since nu is undefined at ijk1,
554     ! approximate gradient with interior points: ijk2 and i,i-1,k relative
555     ! to ijk2
556                 IJK3 = WEST_OF(IJK2)
557                 GNUWDOTN(MM) = 1.d0*(6.d0/(PI*DP_avg(MM)))*&
558                      (EP_s(IJK2,MM)-EP_s(IJK3,MM))*oDX_E(I_OF(IJK3))
559     
560     ! gradient in granular temperature at wall (i+1/2,j,k relative to ijk2)
561     ! dot with the normal vector at the wall.
562                 GTWDOTN(MM) = 1.d0*(Theta_m(IJK2,MM)-Theta_m(IJK3,MM))*&
563                      oDX_E(I_OF(IJK3))
564              ENDDO
565     
566     ! Calculate velocity components at i+1/2, j, k (relative to IJK2)
567              UGC  = U_g(IJK2)
568              VGC  = AVG_X(AVG_Y_N(V_g(JM_OF(IJK2)), V_g(IJK2)),&
569                          AVG_Y_N(V_g(JM_OF(IJK1)), V_g(IJK1)),&
570                          I_OF(IJK2))
571              WGC  = AVG_X(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
572                          AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
573                          I_OF(IJK2))
574              USCM = U_s(IJK2,M)
575              VSCM = AVG_X(AVG_Y_N(V_s(JM_OF(IJK2),M), V_s(IJK2,M)),&
576                          AVG_Y_N(V_s(JM_OF(IJK1),M), V_s(IJK1,M)),&
577                          I_OF(IJK2))
578              WSCM = AVG_X(AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
579                          AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
580                          I_OF(IJK2))
581     
582     ! slip velocity at the wall
583              VSLIPSQ=(WSCM-BC_Ww_s(L, M))**2 + (VSCM-BC_Vw_s(L, M))**2
584     
585     
586           ELSEIF(FCELL== 'T')THEN
587              DO MM = 1, SMAX
588                 TH_avg(MM) = AVG_Z(Theta_m(IJK1,MM),Theta_m(IJK2,MM),K_OF(IJK1))
589                 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta ! for some corner cells
590     
591     ! added for IA (2005) theory:
592     ! include -1 since normal vector is pointing in bottom dir (-z)
593     ! velocity at wall (i,j,k+1/2 relative to ijk1) dot with the normal
594     ! vector at the wall.
595                 VWDOTN(MM) = -1.d0*W_S(IJK1,MM)
596     
597     ! gradient in number density at wall (i,j,k+1/2 relative to ijk1) dot
598     ! with the normal vector at the wall. since nu is undefined at ijk1,
599     ! approximate gradient with interior points: ijk2 and i,j,k+1 relative
600     ! to ijk2
601                 IJK3 = TOP_OF(IJK2)
602                 GNUWDOTN(MM) = -1.d0*(6.d0/(PI*DP_avg(MM)))*&
603                      (EP_s(IJK3,MM)-EP_s(IJK2,MM))*oX(I_of(IJK2))*oDZ_T(K_OF(IJK2))
604     
605     ! gradient in granular temperature at wall (i,j,k+1/2 relative to ijk1)
606     ! dot with the normal vector at the wall.
607                 GNUWDOTN(MM) = -1.d0*(Theta_m(IJK3,MM)-Theta_m(IJK2,MM))*&
608                      oX(I_of(IJK2))*oDZ_T(K_OF(IJK2))
609              ENDDO
610     
611     ! Calculate velocity components at i, j, k+1/2 (relative to IJK1)
612              UGC  = AVG_Z(AVG_X_E(U_g(IM_OF(IJK1)),U_g(IJK1),I_OF(IJK1)),&
613                          AVG_X_E(U_g(IM_OF(IJK2)),U_g(IJK2),I_OF(IJK2)),&
614                          K_OF(IJK1))
615              VGC  = AVG_Z(AVG_Y_N(V_g(JM_OF(IJK1)), V_g(IJK1)),&
616                          AVG_Y_N(V_g(JM_OF(IJK2)), V_g(IJK2)),&
617                          K_OF(IJK1))
618              WGC  = W_g(IJK1)
619              USCM = AVG_Z(AVG_X_E(U_s(IM_OF(IJK1),M),U_s(IJK1,M),I_OF(IJK1)),&
620                          AVG_X_E(U_s(IM_OF(IJK2),M),U_s(IJK2,M),I_OF(IJK2)),&
621                          K_OF(IJK1))
622              VSCM = AVG_Z(AVG_Y_N(V_s(JM_OF(IJK1),M), V_s(IJK1,M)),&
623                          AVG_Y_N(V_s(JM_OF(IJK2),M), V_s(IJK2,M)),&
624                          K_OF(IJK1))
625              WSCM = W_s(IJK1,M)
626     
627     ! slip velocity at the wall
628              VSLIPSQ=(VSCM-BC_Vw_s(L, M))**2 + (USCM-BC_Uw_s(L, M))**2
629     
630     
631           ELSEIF(FCELL== 'B')THEN
632              DO MM = 1, SMAX
633                 TH_avg(MM) = AVG_Z(Theta_m(IJK2,MM),Theta_m(IJK1,MM),K_OF(IJK2))
634                 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta ! for some corner cells
635     
636     ! added for IA (2005) theory:
637     ! include 1 since normal vector is pointing in bottom dir (+z)
638     ! velocity at wall (i,j,k+1/2 relative to ijk2) dot with the normal
639     ! vector at the wall
640                 VWDOTN(MM) = 1.d0*W_S(IJK2,MM)
641     
642     ! gradient in number density at wall (i,j,k+1/2 relative to ijk2) dot
643     ! with the normal vector at the wall. since nu is undefined at ijk1,
644     ! approximate gradient with interior points: ijk2 and i,j,k-1 relative
645     ! to ijk2
646                 IJK3 = BOTTOM_OF(IJK2)
647                 GNUWDOTN(MM) = 1.d0*(6.d0/(PI*DP_avg(MM)))*&
648                      (EP_s(IJK2,MM)-EP_s(IJK3,MM))*oX(I_of(IJK3))*oDZ_T(K_OF(IJK3))
649     
650     ! gradient in granular temperature at wall (i,j,k+1/2 relative to ijk2)
651     ! dot with the normal vector at the wall.
652                 GTWDOTN(MM) = 1.d0*(Theta_m(IJK2,MM)-Theta_m(IJK3,MM))*&
653                      oX(I_of(IJK3))*oDZ_T(K_OF(IJK3))
654              ENDDO
655     
656     ! Calculate velocity components at i, j, k+1/2 (relative to IJK2)
657              UGC  = AVG_Z(AVG_X_E(U_g(IM_OF(IJK2)),U_g(IJK2),I_OF(IJK2)),&
658                          AVG_X_E(U_g(IM_OF(IJK1)),U_g(IJK1),I_OF(IJK1)),&
659                          K_OF(IJK2))
660              VGC  = AVG_Z(AVG_Y_N(V_g(JM_OF(IJK2)), V_g(IJK2)),&
661                          AVG_Y_N(V_g(JM_OF(IJK1)), V_g(IJK1)),&
662                          K_OF(IJK2))
663              WGC  = W_g(IJK2)
664              USCM = AVG_Z(AVG_X_E(U_s(IM_OF(IJK2),M),U_s(IJK2,M),I_OF(IJK2)),&
665                          AVG_X_E(U_s(IM_OF(IJK1),M),U_s(IJK1,M),I_OF(IJK1)),&
666                          K_OF(IJK2))
667              VSCM = AVG_Z(AVG_Y_N(V_s(JM_OF(IJK2),M), V_s(IJK2,M)),&
668                          AVG_Y_N(V_s(JM_OF(IJK1),M), V_s(IJK1,M)),&
669                          K_OF(IJK2))
670              WSCM = W_s(IJK2,M)
671     
672     ! slip velocity at the wall
673              VSLIPSQ=(VSCM-BC_Vw_s(L, M))**2 + (USCM-BC_Uw_s(L, M))**2
674     
675           ELSE
676              WRITE(LINE,'(A, A)') 'Error: Unknown FCELL'
677              CALL WRITE_ERROR('CALC_THETA_BC', LINE, 1)
678              call exitMPI(myPE)
679           ENDIF
680     
681     ! magnitude of gas-solids relative velocity
682           VREL = DSQRT( (UGC-USCM)**2 + (VGC-VSCM)**2 + (WGC-WSCM)**2 )
683     
684           CALL THETA_Hw_Cw(g0, EPs_avg, EPg_avg, ep_star_avg, &
685                            g0EPs_avg, TH_avg,Mu_g_avg,RO_g_avg, ROs_avg, &
686                            DP_avg, K_12_avg,Tau_12_avg,VREL,VSLIPSQ,VWDOTN,&
687                            GNUWDOTN,GTWDOTN,M,Gw,Hw,Cw,L)
688     
689     ! Output the variable specularity coefficient.  Currently only works
690     ! for one solids phase
691           IF(BC_JJ_M .and. PHIP_OUT_JJ) THEN
692              IF(PHIP_OUT_ITER.eq.1) THEN
693                 ReactionRates(IJK1,1)= PHIP_JJ(dsqrt(VSLIPSQ),TH_avg(1))
694              ENDIF
695           ENDIF
696     
697           RETURN
698           END SUBROUTINE CALC_THETA_BC
699     
700     
701     
702     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
703     !                                                                      C
704     !  Subroutine: SUBROUTINE THETA_HW_CW                                  C
705     !  Purpose: Subroutine for gw, hw and cw                               C
706     !                                                                      C
707     !  Author: Kapil Agrawal, Princeton University         Date: 15-MAR-98 C
708     !  Reviewer:                                           Date:           C
709     !                                                                      C
710     !                                                                      C
711     !  Modified: Sofiane Benyahia, Fluent Inc.             Date: 02-FEB-05 C
712     !  Purpose: Include conductivity defined by Simonin and Ahmadi         C
713     !           Also included Jenkins small frictional limit               C
714     !                                                                      C
715     !  Literature/Document References:                                     C
716     !     See calc_mu_s.f for references on kinetic theory models          C
717     !     Johnson, P. C., and Jackson, R., Frictional-collisional          C
718     !        constitutive relations for granular materials, with           C
719     !        application to plane shearing, Journal of Fluid Mechanics,    C
720     !        1987, 176, pp. 67-93                                          C
721     !     Jenkins, J. T., and Louge, M. Y., On the flux of fluctuating     C
722     !        energy in a collisional grain flow at a flat frictional wall, C
723     !        Physics of Fluids, 1997, 9(10), pp. 2835-2840                 C
724     !                                                                      C
725     !                                                                      C
726     !  Additional Notes:                                                   C
727     !     The JJ granular energy BC is written as: n.q = D-G, where        C
728     !       q = the heat flux, n = outward normal vector                   C
729     !       D = dissipation due to p-w collisions, and                     C
730     !       G = generation due to p-w slip                                 C
731     !     In MFIX this is implemented as k.gradT = G-D, where              C
732     !       k = granular conductivity units (Mass/Length/Time)             C
733     !       T = granular temperature in units (Length/Time)^2              C
734     !     However, the expression for heat flux (q) may contain terms      C
735     !     beyond that of gradient in temperature, such as the dufour       C
736     !     term. A more rigorous implemenation should account for these     C
737     !     additional terms, which is not done here.                        C
738     !                                                                      C
739     !                                                                      C
740     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
741     
742           SUBROUTINE THETA_HW_CW(g0,EPS,EPG, ep_star_avg, &
743                                  g0EPs_avg,TH,Mu_g_avg,RO_g_avg, ROs_avg, DP_avg, &
744                                  K_12_avg,Tau_12_avg,VREL,VSLIPSQ,VWDOTN,&
745                                  GNUWDOTN,GTWDOTN,M,GW,HW,CW,L)
746     
747     !-----------------------------------------------
748     ! Modules
749     !-----------------------------------------------
750           USE bc
751           USE compar
752           USE constant
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     
1284     
1285              CASE DEFAULT
1286     ! should never hit this
1287                 WRITE (*, '(A)') 'BC_THETA => THETA_HW_CW'
1288                 WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1289                 call mfix_exit(myPE)
1290              END SELECT
1291     
1292     
1293           RETURN
1294           END SUBROUTINE THETA_HW_CW
1295     
1296     
1297     
1298     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1299     !                                                                      C
1300     !  Function: PHIP_JJ                                                   C
1301     !  Purpose: Calculate the specularity coefficient for JJ BC            C
1302     !                                                                      C
1303     !  Author: Tingwen Li                                                  C
1304     !                                                                      C
1305     !                                                                      C
1306     !  Literature/Document References:                                     C
1307     !     Li, T., and Benyahia, S., Revisiting Johnson and Jackson         C
1308     !        boundary conditions for granular flows, AIChE Journal, 2012,  C
1309     !        Vol 58 (7), pp. 2058-2068                                     C
1310     !                                                                      C
1311     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1312     
1313     
1314           DOUBLE PRECISION FUNCTION phip_jj(uslip,g_theta)
1315     !-----------------------------------------------
1316     ! Modules
1317     !-----------------------------------------------
1318           use constant, only : PI, k4phi, phip0
1319           implicit none
1320     !-----------------------------------------------
1321     ! Dummy arguments
1322     !-----------------------------------------------
1323           double precision, INTENT(IN) :: uslip
1324           double precision, INTENT(IN) :: g_theta
1325     !-----------------------------------------------
1326     ! Local variables
1327     !-----------------------------------------------
1328           double precision :: r4phi
1329     !-----------------------------------------------
1330     
1331           !k4phi=7.d0/2.d0*mu4phi*(1.d0+e_w)
1332           r4phi=uslip/dsqrt(3.d0*g_theta)
1333     
1334           IF(r4phi .le. 4.d0*k4phi/7.d0/dsqrt(6.d0*PI)/phip0) THEN
1335              phip_jj=-7.d0*dsqrt(6.d0*PI)*phip0**2/8.d0/k4phi*r4phi+phip0
1336           ELSE
1337              phip_jj=2.d0/7.d0*k4phi/r4phi/dsqrt(6.d0*PI)
1338           ENDIF
1339     
1340           RETURN
1341     
1342           END FUNCTION PHIP_JJ
1343