File: N:\mfix\model\calc_mu_s.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_MU_s                                               C
4     !  Purpose: Calculate the vicosity, second viscosity, and pressure     C
5     !  associated with the Mth 'solids' phase. In the 'default' or         C
6     !  undefined case closure for 'granular conductivity' is also          C
7     !  calculated.                                                         C
8     !                                                                      C
9     !                                                                      C
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
11           SUBROUTINE CALC_MU_s(M)
12     
13     ! Modules
14     !---------------------------------------------------------------------//
15     ! some constants
16           USE param1, only: undefined
17     ! specified constant solids phase viscosity
18           USE physprop, only: mu_s0
19     ! invoke user defined quantity
20           USE usr_prop, only: usr_mus, calc_usr_prop
21           USE usr_prop, only: solids_viscosity
22           IMPLICIT NONE
23     
24     ! Dummy arguments
25     !---------------------------------------------------------------------//
26     ! solids phase index
27           INTEGER, INTENT(IN) :: M
28     
29     !---------------------------------------------------------------------//
30     
31           IF (USR_MUS(M)) THEN
32              CALL CALC_USR_PROP(Solids_Viscosity,lm=M)
33           ELSEIF (MU_S0(M) == UNDEFINED) THEN
34     ! this is a necessary check as one may have mu_s0 defined but still 
35     ! need to call this routine (for set_epmus)
36              CALL CALC_DEFAULT_MUs(M)
37           ENDIF
38     
39           CALL SET_EPMUS_VALUES(M)
40           RETURN
41           END SUBROUTINE CALC_MU_s
42     
43     
44     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
45     !                                                                      C
46     !  Subroutine: SET_EPMUS_VALUES                                        C
47     !  Purpose: This routine sets the internal variables epmu_s and        C
48     !  eplambda_s that are used in the stress calculations. If the         C
49     !  keyword Ishii is invoked then these quantities represent the        C
50     !  viscosity and second viscosity multiplied by the volume fraction    C
51     !  otherwise they are simply viscosity/second viscosity (i.e. are      C
52     !  multipled by one).                                                  C
53     !                                                                      C
54     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
55           SUBROUTINE SET_EPMUS_VALUES(M)
56     
57     ! Modules
58     !---------------------------------------------------------------------//
59           USE compar, only: ijkstart3, ijkend3
60           USE fldvar, only: eps_ifac
61           USE functions, only: fluid_at
62           USE param1, only: zero
63           USE visc_s, only: mu_s, epmu_s, lambda_s, eplambda_s
64           USE mms, only: use_mms
65     
66     ! Dummy arguments
67     !---------------------------------------------------------------------//
68     ! solids phase index
69           INTEGER, INTENT(IN) :: M
70     
71     ! Local variables
72     !---------------------------------------------------------------------//
73     ! cell index
74           INTEGER :: IJK
75     !---------------------------------------------------------------------//
76     
77     !      EPMU_S(:,M) = EPS_IFAC(:,M)*MU_s(:,M)
78     !      EPLAMBDA_S(:,M) = EPS_IFAC(:,M)*LAMBDA_S(:,M)
79     
80           DO IJK=ijkstart3,ijkend3
81              IF(FLUID_AT(IJK) .OR. USE_MMS) THEN
82     ! if ishii then multiply by volume fraction otherwise multiply by 1
83                 EPMU_S(IJK,M) = EPS_IFAC(IJK,M)*Mu_s(IJK,M)
84                 EPLAMBDA_S(IJK,M) = EPS_IFAC(IJK,M)*Lambda_s(IJK,M)
85              ELSE
86                 EPMU_S(IJK,M) = ZERO
87                 EPLAMBDA_S(IJK,M) = ZERO
88              ENDIF
89           ENDDO
90           RETURN
91           END SUBROUTINE SET_EPMUS_VALUES
92     
93     
94     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
95     !                                                                      C
96     !  Subroutine CALC_DEFAULT_MUS                                         C
97     !  Purpose: By default MFIX assumes each Mth phase is a granular       C
98     !  solids material. Accordingly, various theories are invoked to       C
99     !  provide closures for the solids stress term. These include          C
100     !  different kinetic theory and soil mechanic models, which may be     C
101     !  specified by the user through various keyword selections.           C
102     !  The granular stress terms (granular viscosity, second/bulk          C
103     !  viscosity, solids pressure), and if required, a granular            C
104     !  conductivity term are closed here.                                  C
105     !                                                                      C
106     !  Comments:                                                           C
107     !     GRANULAR_ENERGY = .FALSE.                                        C
108     !        EP_g < EP_star   -->    friction_schaeffer                    C
109     !        EP_g >= EP_star  -->    viscous (algebraic)                   C
110     !                                                                      C
111     !     GRANULAR_ENERGY = .TRUE.                                         C
112     !        FRICTION = .TRUE.                                             C
113     !           EP_s(IJK,M) > EPS_f_min  -->  friction + viscous(pde)      C
114     !           EP_s(IJK,M) < EP_f_min   -->  viscous (pde)                C
115     !        FRICTION = .FALSE.                                            C
116     !           EP_g < EP_star  -->  friction_schaeffer + viscous(pde)     C
117     !           EP_g >= EP_star -->  viscous (pde)                         C
118     !                                                                      C
119     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
120           SUBROUTINE CALC_DEFAULT_MUS(M)
121     
122     ! Modules
123     !---------------------------------------------------------------------//
124           USE compar, only: ijkstart3, ijkend3
125     !
126           USE error_manager, only: err_msg, init_err_msg, finl_err_msg
127           USE error_manager, only: flush_err_msg
128     ! solids pressure
129           USE fldvar, only: p_s, p_s_c, p_s_v
130           USE fldvar, only: p_s_f
131     ! some constants
132           USE param1, only: one
133     ! runtime flag indicates whether the solids phase becomes close-packed
134     ! at ep_star
135           USE physprop, only: close_packed
136     ! number of solids phases
137           USE physprop, only: mmax
138           USE physprop, only: blend_function
139     
140     ! runtime flag to use qmomk
141           USE qmom_kinetic_equation, only: qmomk
142     
143     ! runtime flag to solve partial differential granular energy eqn(s)
144           USE run, only: granular_energy
145     ! kinetic theories
146           USE run, only: kt_type_enum
147           USE run, only: lun_1984
148           USE run, only: simonin_1996
149           USE run, only: ahmadi_1995
150           USE run, only: gd_1999
151           USE run, only: gtsh_2012
152           USE run, only: ia_2005
153           USE run, only: ghd_2007
154           USE run, only: kt_type
155     ! filtered subgrid model
156           USE run, only: subgrid_type_enum
157           USE run, only: milioli, igci, undefined_subgrid_type
158     ! frictional theories
159           USE run, only: friction, schaeffer
160     ! runtime flag for blending stress (only with schaeffer)
161           USE run, only: blending_stress
162     ! solids transport coefficients
163           USE visc_s, only: mu_s, epmu_s, mu_s_c, mu_s_v
164           USE visc_s, only: mu_s_p, mu_s_f
165           USE visc_s, only: lambda_s, eplambda_s, lambda_s_c, lambda_s_v
166           USE visc_s, only: lambda_s_p, lambda_s_f
167           IMPLICIT NONE
168     
169     ! Dummy arguments
170     !---------------------------------------------------------------------//
171     ! solids phase index
172           INTEGER, INTENT(IN) :: M
173     
174     ! Local variables
175     !---------------------------------------------------------------------//
176     ! cell index
177           INTEGER :: IJK
178     ! blend factor
179           DOUBLE PRECISION :: BLEND
180     
181     ! General initializations
182     !---------------------------------------------------------------------
183     ! Zero out various quantities
184           CALL INIT0_MU_S(M)
185     
186     ! Calculate quantities that are functions of velocity only and may be
187     ! needed for subsequent calculations by various model options
188           CALL INIT1_MU_S(M)
189     
190     
191     ! Viscous-flow stress tensor
192     !---------------------------------------------------------------------
193           IF (.NOT. QMOMK) THEN
194     ! if QMOMK then do not solve algebraic or PDE form of granular
195     ! temperature governing equation
196              IF(.NOT.GRANULAR_ENERGY) THEN
197                 IF(SUBGRID_TYPE_ENUM .ne. UNDEFINED_SUBGRID_TYPE) THEN
198                    IF (SUBGRID_TYPE_ENUM .EQ. IGCI) THEN
199                       CALL SUBGRID_STRESS_IGCI(M)
200                    ELSEIF (SUBGRID_TYPE_ENUM .EQ. MILIOLI) THEN
201                       CALL SUBGRID_STRESS_MILIOLI(M)
202                    ENDIF
203                 ELSE
204                   call gt_algebraic(M)   ! algebraic granular energy equation
205                 ENDIF
206              ELSE   ! granular energy transport equation
207                 SELECT CASE(KT_TYPE_ENUM)
208                    CASE (IA_2005)   ! polydisperse theory
209                       CALL gt_pde_ia(M)
210                    CASE (GD_1999)   ! strictly monodisperse theory
211                       CALL gt_pde_gd(M)
212                    CASE (GTSH_2012)   ! strictly monodisperse theory
213                       CALL gt_pde_gtsh(M)
214                    CASE (GHD_2007)   ! polydisperse GHD theory for mixture temperature
215                       CALL TRANSPORT_COEFF_GHD(M)
216                    CASE(LUN_1984)   ! monodisperse/ad-hoc polydisperse theory
217                       CALL gt_pde_lun(M)
218                    CASE(SIMONIN_1996)   ! monodisperse/ad-hoc polydisperse theory
219                       CALL gt_pde_simonin(M)
220                    CASE (AHMADI_1995)   ! monodisperse/ad-hoc polydisperse theory
221                       CALL gt_pde_ahmadi(M)
222                    CASE DEFAULT
223     ! should never hit this
224                       CALL INIT_ERR_MSG("CALC_MU_S")
225                       WRITE(ERR_MSG, 1100) KT_TYPE
226      1100 FORMAT('ERROR 1100: Invalid KT_TYPE: ', A,' The check_data ',&
227              'routines should',/, 'have already caught this error and ',&
228              'prevented the simulation from ',/,'running. Please notify ',&
229              'the MFIX developers.')
230                       CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
231                       CALL FINL_ERR_MSG
232                 END SELECT  ! end selection of kt_type_enum
233              ENDIF
234           ENDIF
235     
236     
237     ! Frictional stress tensors
238     ! Note that only one of these can be used at this time
239     !---------------------------------------------------------------------
240     ! Schaeffer's frictional formulation
241           IF (SCHAEFFER .AND. CLOSE_PACKED(M)) call friction_schaeffer(M)
242     ! Princeton's frictional implementation
243           IF (FRICTION .AND. CLOSE_PACKED(M)) call friction_princeton(M)
244     
245     
246     ! Assign viscosity, second viscosity and pressure in mth solids phase.
247     ! Note that a plastic pressure component is calculated in a separate
248     ! routine (see calc_p_star) which is then directly incorporated into
249     ! the solids momentum equations (see source_u_s, source_v_s and
250     ! source_w_s).
251     !---------------------------------------------------------------------
252           Mu_s_c(:,M) = Mu_s_v(:)
253           LAMBDA_s_c(:,M)= Lambda_s_v(:)
254           P_s_c(:,M) = P_s_v(:)
255     
256           IF(BLENDING_STRESS) THEN
257     ! blend plastic & viscous stresses (not available with friction)
258              DO IJK = ijkstart3, ijkend3
259                 blend =  blend_function(IJK)
260                 Mu_s(IJK,M) = (ONE-blend)*Mu_s_p(IJK) &
261                     + blend*Mu_s_v(IJK)
262     ! is there any point in blending P_s_p or lambda_s_p since neither are
263     ! assigned? the only thing this would effectively do is reduce P_s_v
264     ! and lambda_s_v by the blend value and therefore P_s and lambda_s.
265                 LAMBDA_s(IJK,M) = (ONE-blend)*Lambda_s_p(IJK) + &
266                    blend*Lambda_s_v(IJK)
267              ENDDO
268           ELSE
269              Mu_s(:,M) = Mu_s_p(:) + Mu_s_v(:) + Mu_s_f(:)
270              LAMBDA_s(:,M) = Lambda_s_p(:) + Lambda_s_v(:) + Lambda_s_f(:)
271              P_s(:,M) = P_s_v(:) + P_s_f(:)
272           ENDIF  ! end if/else (blending_stress)
273     
274           RETURN
275           END SUBROUTINE CALC_DEFAULT_MUs
276     
277     
278     
279     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
280     !                                                                      C
281     !  Subroutine: gt_algebraic                                            C
282     !  Purpose: solve an algebraic form of the granular energy equation    C
283     !  based on Lun et al. (1984)                                          C
284     !  Obtained from the granular energy equation of Lun et al. (1984) by  C
285     !  assuming that the granular energy is dissipated locally so that     C
286     !  diffusion and convection contributions are neglected and only       C
287     !  generation and dissipation terms are retained.                      C
288     !                                                                      C
289     !  References:                                                         C
290     !  Syamlal, M., 1987c, "A Review of Granular Stress Constitutive       C
291     !     Relations," Topical Report, DOE/MC/21353-2372, NTIS/DE87006499,  C
292     !     National Technical Information Service, Springfield, VA.         C
293     !                                                                      C
294     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
295           Subroutine Gt_algebraic (M)
296     
297     ! Modules
298     !---------------------------------------------------------------------//
299           USE compar, only: ijkstart3, ijkend3
300           USE constant, only: C_e, sqrt_pi
301           USE constant, only: v_ex
302           USE fldvar, only: p_s_v
303           USE fldvar, only: ep_g, ep_s
304           USE fldvar, only: d_p, ro_s
305           USE fldvar, only: theta_m
306           USE functions, only: fluid_at
307           USE param1, only: zero, half, one, small_number
308           USE rdf, only: g_0
309           USE trace, only: trd_s_c, trd_s2
310           USE visc_s, only: mu_s_v, lambda_s_v
311           use visc_s, only: ep_g_blend_start
312           USE visc_s, only: alpha_s
313           USE visc_s, only: trm_s, trdM_s
314           IMPLICIT NONE
315     
316     ! Dummy arguments
317     !---------------------------------------------------------------------//
318     ! solids phase index
319           INTEGER, INTENT(IN) :: M
320     
321     ! Local variables
322     !---------------------------------------------------------------------//
323     ! cell index
324           INTEGER :: IJK
325     ! Coefficients of quadratic equation
326           DOUBLE PRECISION :: aq, bq, cq
327     ! Constant in equation for mth solids phase pressure
328           DOUBLE PRECISION :: K_1m
329     ! Constant in equation for mth solids phase bulk viscosity
330           DOUBLE PRECISION :: K_2m
331     ! Constant in equation for mth solids phase viscosity
332           DOUBLE PRECISION :: K_3m
333     ! Constant in equation for mth solids phase dissipation
334           DOUBLE PRECISION :: K_4m
335     ! Constant in Boyle-Massoudi stress term
336           DOUBLE PRECISION :: K_5m
337     ! Value of EP_s * SQRT( THETA )for Mth solids phase continuum
338           DOUBLE PRECISION :: EP_sxSQRTHETA
339     ! Value of EP_s * EP_s * THETA for Mth solids phase continuum
340           DOUBLE PRECISION :: EP_s2xTHETA, temp_local
341     !---------------------------------------------------------------------//
342     
343     !!$omp parallel do default(shared)                                    &
344     !!$omp private(IJK, K_1m, K_2m, K_3m, K_4m, K_5m, temp_local,         &
345     !!$omp         aq, bq, cq, EP_sxSQRTHETA, EP_s2xTHETA)
346            DO IJK = ijkstart3, ijkend3
347              IF ( FLUID_AT(IJK) ) THEN
348                 IF(EP_g(IJK) .GE. EP_g_blend_start(IJK)) THEN
349     
350     ! Calculate K_1m, K_2m, K_3m, K_4m
351                    K_1m = 2.D0 * (ONE + C_e) * RO_S(IJK,M) * G_0(IJK, M, M)
352                    K_3m = HALF * D_p(IJK,M) * RO_S(IJK,M) * ( ( (SQRT_PI / &
353                       (3.D0 * (3.D0 - C_e))) * ( HALF*(3d0*C_e+ONE) + &
354                       0.4D0*(ONE + C_e)*(3.D0*C_e - ONE) * EP_s(IJK,M)* &
355                       G_0(IJK, M,M)) ) + 8.D0*EP_s(IJK,M)*G_0(IJK, M,M)*&
356                       (ONE + C_e)/ (5.D0*SQRT_PI) )
357                    K_2m = 4.D0 * D_p(IJK,M) * RO_S(IJK,M) * (ONE + C_e) *&
358                       EP_s(IJK,M) * G_0(IJK, M,M) / (3.D0 * SQRT_PI) - &
359                       2.D0/3.D0 * K_3m
360                    K_4m = 12.D0 * (ONE - C_e*C_e) *&
361                       RO_S(IJK,M) * G_0(IJK, M,M) / (D_p(IJK,M) * SQRT_PI)
362                    aq = K_4m*EP_s(IJK,M)
363                    bq = K_1m*EP_s(IJK,M)*trD_s_C(IJK,M)
364                    cq = -(K_2m*trD_s_C(IJK,M)*trD_s_C(IJK,M)&
365                           + 2.D0*K_3m*trD_s2(IJK,M))
366     
367     ! Boyle-Massoudi Stress term
368                    IF(V_ex .NE. ZERO) THEN
369                       K_5m = 0.4 * (ONE + C_e) * G_0(IJK, M,M) * &
370                          RO_S(IJK,M) * ( (V_ex * D_p(IJK,M)) / &
371                          (ONE - EP_s(IJK,M) * V_ex) )**2
372                       bq = bq + EP_s(IJK,M) * K_5m * &
373                          (trM_s(IJK) + 2.D0 * trDM_s(IJK))
374                    ELSE
375                       K_5m = ZERO
376                    ENDIF
377     
378     ! Calculate EP_sxSQRTHETA and EP_s2xTHETA
379                    temp_local = bq**2 - 4.D0 * aq * cq
380                    EP_sxSQRTHETA = (-bq + SQRT(temp_local))&
381                       / ( 2.D0 * K_4m )
382                    EP_s2xTHETA = EP_sxSQRTHETA * EP_sxSQRTHETA
383     
384                    IF(EP_s(IJK,M) > SMALL_NUMBER)THEN
385     ! Find pseudo-thermal temperature in the Mth solids phase
386                       THETA_m(IJK,M) = EP_s2xTHETA/(EP_s(IJK,M)*EP_s(IJK,M))
387                    ELSE
388                       THETA_m(IJK,M) = ZERO
389                    ENDIF
390     
391     ! Find pressure in the Mth solids phase
392                    P_s_v(IJK) = K_1m * EP_s2xTHETA
393     
394     ! second viscosity in Mth solids phase
395                    LAMBDA_s_v(IJK) = K_2m * EP_sxSQRTHETA
396     
397     ! shear viscosity in Mth solids phase
398                    MU_s_v(IJK) = K_3m * EP_sxSQRTHETA
399     
400     ! Boyle-Massoudi stress coefficient
401                    ALPHA_s(IJK, M) = -K_5m * EP_s2xTHETA
402                 ENDIF
403     
404              ENDIF   ! Fluid_at
405           ENDDO
406     !!$omp end parallel do
407     
408           RETURN
409           END SUBROUTINE GT_ALGEBRAIC
410     
411     
412     
413     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
414     !                                                                      C
415     !  Subroutine: GT_PDE_LUN                                              C
416     !  Purpose: Calculate granular stress terms (viscosity, bulk viscosity C
417     !  solids pressure) & granular conductivity                            C
418     !                                                                      C
419     !  Author: Kapil Agrawal, Princeton University      Date: 6-FEB-98     C
420     !                                                                      C
421     !  Literature/Document References:                                     C
422     !  Lun, C.K.K., S.B. Savage, D.J. Jeffrey, and N. Chepurniy,           C
423     !     Kinetic theories for granular flow - inelastic particles in      C
424     !     Couette-flow and slightly inelastic particles in a general       C
425     !     flow field. Journal of Fluid Mechanics, 1984. 140(MAR):          C
426     !     p. 223-256                                                       C
427     !                                                                      C
428     !                                                                      C
429     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
430           Subroutine gt_pde_lun (M)
431     
432     ! Modules
433     !---------------------------------------------------------------------//
434           USE constant, only: switch
435           USE constant, only: alpha
436           USE constant, only: pi
437           USE constant, only: eta
438     
439           USE drag, only: f_gs
440     
441           USE fldvar, only: ro_g
442           USE fldvar, only: ep_s
443           USE fldvar, only: d_p, rop_s, ro_s
444           USE fldvar, only: theta_m
445           USE fldvar, only: p_s_v
446     
447           USE kintheory, only: kt_dga
448     
449           USE param1, only: zero, one, small_number
450     
451           USE physprop, only: smax
452           USE physprop, only: Kth_s, kphi_s
453     
454           USE rdf, only: g_0, DG_0DNU
455     
456           USE toleranc, only: dil_ep_s
457     
458           USE visc_s, only: mu_s_v, mu_b_v, lambda_s_v
459     
460           USE compar, only: ijkstart3, ijkend3
461           USE functions, only: fluid_at
462           IMPLICIT NONE
463     
464     ! Dummy arguments
465     !---------------------------------------------------------------------//
466     ! solids phase index
467           INTEGER, INTENT(IN) :: M
468     
469     ! Local variables
470     !---------------------------------------------------------------------//
471     ! cell indices
472           INTEGER :: IJK
473     ! solids phase index
474           INTEGER :: MM
475     ! use to compute MU_s(IJK,M) & Kth_S(IJK,M)
476           DOUBLE PRECISION :: Mu_star, Kth, Kth_star
477     ! sum of ep_s * g_0
478           DOUBLE PRECISION :: SUM_EpsGo
479     ! single particle drag coefficient/ep_s
480           DOUBLE PRECISION :: dga_sm
481     !---------------------------------------------------------------------//
482     
483           DO IJK = ijkstart3, ijkend3
484              IF ( FLUID_AT(IJK) ) THEN
485     
486     ! The following is purely an ad-hoc modification so that the underlying
487     ! monodisperse theory can be used for polydisperse systems in a
488     ! consistent manner.  That is, solids pressure, viscosity and
489     ! conductivity must be additive. The non-linear terms (eps^2) are
490     ! corrected so the stresses of two or more identical solids phases are
491     ! equal to those of a equivalent single solids phase. sof June 15 2005.
492                 SUM_EpsGo = ZERO
493                 DO MM = 1, SMAX
494                    SUM_EpsGo =  SUM_EpsGo+EP_s(IJK,MM)*G_0(IJK,M,MM)
495                 ENDDO
496     
497                 IF(EP_S(IJK,M) < DIL_EP_S) &
498                    dga_sM = KT_DGA(ijk, M)
499     
500     ! Pressure
501                 P_s_v(IJK) = ROP_s(IJK,M)*(1.d0 + 4.d0 * Eta *&
502                     SUM_EpsGo)*Theta_m(IJK,M)
503     
504     ! Bulk and shear viscosity
505                 Mu_s_v(IJK) = (5.d0*DSQRT(Pi*Theta_m(IJK,M))*D_p(IJK,M)*&
506                    RO_S(IJK,M))/96.d0
507                 Mu_b_v(IJK) = (256.d0*Mu_s_v(IJK)*EP_s(IJK,M)*SUM_EpsGo)&
508                     /(5.d0*Pi)
509     
510     ! added Ro_g = 0 for granular flows (no gas). sof Aug-02-2005
511                 IF(SWITCH == ZERO .OR. RO_G(IJK) == ZERO) THEN
512                    Mu_star = Mu_s_v(IJK)
513                 ELSEIF(Theta_m(IJK,M) .LT. SMALL_NUMBER)THEN
514                    Mu_star = ZERO
515                 ELSEIF(EP_S(IJK,M) < DIL_EP_S) THEN
516                    Mu_star = RO_S(IJK,M)*EP_s(IJK,M)* G_0(IJK,M,M)*&
517                        Theta_m(IJK,M)* Mu_s_v(IJK)/ &
518                        (RO_S(IJK,M)*SUM_EpsGo*Theta_m(IJK,M) &
519                        + 2.d0*DgA_SM/RO_S(IJK,M)* Mu_s_v(IJK))
520                 ELSE
521                    Mu_star = RO_S(IJK,M)*EP_s(IJK,M)* G_0(IJK,M,M)*&
522                        Theta_m(IJK,M)*Mu_s_v(IJK)/ &
523                        (RO_S(IJK,M)*SUM_EpsGo*Theta_m(IJK,M)+ &
524                        (2.d0*F_gs(IJK,M)*Mu_s_v(IJK)/&
525                        (RO_S(IJK,M)*EP_s(IJK,M) )) )
526                 ENDIF
527     
528     ! Shear viscosity
529                 Mu_s_v(IJK) =&
530                     ((2.d0+ALPHA)/3.d0)*((Mu_star/(Eta*(2.d0-Eta)*&
531                     G_0(IJK,M,M)))*(ONE+1.6d0*Eta*SUM_EpsGo)&
532                     *(ONE+1.6d0*Eta*(3d0*Eta-2d0)*&
533                     SUM_EpsGo)+(0.6d0*Mu_b_v(IJK)*Eta))
534     
535     ! Second viscosity as defined in MFIX
536                 LAMBDA_S_V(IJK) = Eta*Mu_b_v(IJK) - (2d0*Mu_s_v(IJK))/3d0
537     
538                 Kth=75d0*RO_S(IJK,M)*D_p(IJK,M)*DSQRT(Pi*Theta_m(IJK,M))/&
539                     (48d0*Eta*(41d0-33d0*Eta))
540     
541                 IF(SWITCH == ZERO .OR. RO_G(IJK) == ZERO) THEN
542                    Kth_star=Kth
543                 ELSEIF(Theta_m(IJK,M) .LT. SMALL_NUMBER)THEN
544                    Kth_star = ZERO
545                 ELSEIF(EP_S(IJK,M) < DIL_EP_S) THEN
546                    Kth_star = RO_S(IJK,M)*EP_s(IJK,M)* &
547                       G_0(IJK,M,M)*Theta_m(IJK,M)* Kth/ &
548                       (RO_S(IJK,M)*SUM_EpsGo*Theta_m(IJK,M) &
549                       + 1.2d0*DgA_sM/RO_S(IJK,M)* Kth)
550                 ELSE
551                    Kth_star = RO_S(IJK,M)*EP_s(IJK,M)* &
552                       G_0(IJK,M,M)*Theta_m(IJK,M)*Kth/ &
553                       (RO_S(IJK,M)*SUM_EpsGo*Theta_m(IJK,M)+ &
554                       (1.2d0*F_gs(IJK,M)*Kth/(RO_S(IJK,M)*EP_s(IJK,M))) )
555                 ENDIF
556     
557     ! Granular conductivity
558                 Kth_s(IJK,M) = Kth_star/G_0(IJK,M,M)*(&
559                    ( ONE+(12d0/5.d0)*Eta*SUM_EpsGo ) * &
560                    ( ONE+(12d0/5.d0)*Eta*Eta*(4d0*Eta-3d0)* SUM_EpsGo ) + &
561                    (64d0/(25d0*Pi))*(41d0-33d0*Eta)*(Eta*SUM_EpsGo)**2 )
562     
563     
564     
565     ! Dufour coefficient: granular 'conductivity' in the Mth solids phase
566     ! associated with gradient in volume fraction
567     !--------------------------------------------------------------------
568     !     Kphi_s has been set to zero.  To activate the feature uncomment the
569     !     following lines and also the lines in source_granular_energy.
570                 Kphi_s(IJK,M) = ZERO
571     !     &            (Kth_star/(G_0(IJK,M,M)))*(12d0/5.)*Eta*(Eta-1.)*
572     !     &            (2.*Eta-1.)*(1.+(12d0/5.)*Eta*EP_s(IJK,M)*
573     !     &            G_0(IJK,M,M))*(EP_s(IJK,M)*
574     !     &            DG_0DNU(EP_s(IJK,M))
575     !     &            + 2*G_0(IJK,M,M))*Theta_m(IJK,M)
576     !--------------------------------------------------------------------
577     
578              ENDIF   ! Fluid_at
579           ENDDO   ! IJK loop
580           RETURN
581           END SUBROUTINE GT_PDE_LUN
582     
583     
584     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
585     !                                                                      C
586     !  Subroutine: GT_PDE_ahmadi                                           C
587     !  Purpose: Calculate granular stress terms (viscosity, bulk viscosity C
588     !  solids pressure) & granular conductivity using ahmadi model         C
589     !                                                                      C
590     !  Author: Sofiane Benyahia, Fluent Inc.            Date: 02-01-05     C
591     !                                                                      C
592     !  Literature/Document References:                                     C
593     !  Cao, J. and Ahmadi, G., 1995, Gas-particle two-phase turbulent      C
594     !     flow in a vertical duct. Int. J. Multiphase Flow, vol. 21,       C
595     !     No. 6, pp. 1203-1228.                                            C
596     !                                                                      C
597     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
598           Subroutine gt_pde_ahmadi (M)
599     
600     ! Modules
601     !---------------------------------------------------------------------//
602           USE constant, only: C_e, eta
603     
604           USE drag, only: f_gs
605     
606           USE fldvar, only: p_s_v
607           USE fldvar, only: ROP_s, ro_s, d_p, theta_m
608           USE fldvar, only: ep_s
609           USE fldvar, only: k_turb_g, e_turb_g
610           USE kintheory, only: kt_dga
611     
612           USE param1, only: zero, half, one, small_number
613     
614           USE physprop, only: smax
615           USE physprop, only: kth_s
616     
617           USE rdf, only: g_0
618     
619           USE toleranc, only: dil_ep_s
620     
621           USE turb, only: tau_1
622     
623           USE visc_s, only: lambda_s_v, mu_s_v, mu_b_v
624           USE visc_s, only: ep_star_array
625     
626           USE functions, only: fluid_at
627           USE compar, only: ijkstart3, ijkend3
628           IMPLICIT NONE
629     
630     ! Dummy arguments
631     !---------------------------------------------------------------------//
632     ! solids phase index
633           INTEGER, INTENT(IN) :: M
634     
635     ! Local parameters
636     !---------------------------------------------------------------------//
637     ! constant in simonin model
638           DOUBLE PRECISION, PARAMETER :: C_mu = 9.0D-02
639     
640     ! Local variables
641     !---------------------------------------------------------------------//
642     ! cell indices
643           INTEGER :: IJK
644     ! solids phase index
645           INTEGER :: MM, L
646     ! defining parameters for Ahmadi
647           DOUBLE PRECISION :: Tau_12_st
648           DOUBLE PRECISION :: Tmp_Ahmadi_Const
649     ! sum of ep_s * g_0
650           DOUBLE PRECISION :: SUM_EpsGo
651     ! single particle drag coefficient/ep_s
652           DOUBLE PRECISION :: dga_sl
653     !---------------------------------------------------------------------//
654     
655           DO IJK = ijkstart3, ijkend3
656              IF ( FLUID_AT(IJK) ) THEN
657     
658     ! Define time scales and constants
659     
660     ! time scale of turbulent eddies
661                 Tau_1(ijk) = 3.d0/2.d0*C_MU*K_Turb_G(IJK)/&
662                    (E_Turb_G(IJK)+small_number)
663     
664     ! NOTE: Some calculations are based explicitly on solids phase 1!
665     ! Parameters based on L=1: tau_12_st....
666                 L = 1
667     
668     ! Particle relaxation time. For very dilute flows avoid singularity
669     ! by redefining the drag as single particle drag
670                 IF(Ep_s(IJK,M) > DIL_EP_S .AND. &
671                    F_GS(IJK,L) > small_number) THEN
672                    Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,L)
673                 ELSE
674                    dgA_sl = kt_dga(IJK, L)
675                    Tau_12_st = RO_S(IJK,M)/DgA_SL
676                 ENDIF
677     
678     
679     ! The following is purely an ad-hoc modification so that the underlying
680     ! monodisperse theory can be used for polydisperse systems in a
681     ! consistent manner.  That is, solids pressure, viscosity and
682     ! conductivity must be additive. THe non-linear terms (eps^2) are
683     ! corrected so the stresses of two or more identical solids phases are
684     ! equal to those of a equivalent single solids phase. sof June 15 2005.
685                 SUM_EpsGo = ZERO
686                 DO MM = 1, SMAX
687                    SUM_EpsGo =  SUM_EpsGo+EP_s(IJK,MM)*G_0(IJK,M,MM)
688                 ENDDO
689     
690                 P_s_v(IJK) = ROP_s(IJK,M)*Theta_m(IJK,M) * &
691                    ( (ONE + 4.d0*SUM_EpsGo ) + HALF*(ONE - C_e*C_e) )
692     
693                 IF(EP_s(IJK,M) < (ONE-EP_star_array(ijk))) THEN
694                    Tmp_Ahmadi_Const = &
695                       ONE/(ONE+ Tau_1(ijk)/Tau_12_st * &
696                       (ONE-EP_s(IJK,M)/(ONE-EP_star_array(ijk)))**3)
697                 ELSE
698                    Tmp_Ahmadi_Const = ONE
699                 ENDIF
700     
701     ! Shear viscosity
702     ! Note that Ahmadi coefficient 0.0853 in C_mu was replaced by 0.1567
703     ! to include 3/2*sqrt(3/2) because K = 3/2 Theta_m
704                 Mu_s_v(IJK) = Tmp_Ahmadi_Const &
705                     *0.1045d0*(ONE/G_0(IJK,M,M)+3.2d0*EP_s(IJK,M)+12.1824d0*   &
706                     G_0(IJK,M,M)*EP_s(IJK,M)*EP_s(IJK,M))*D_p(IJK,M)*RO_S(IJK,M)*  &
707                     DSQRT(Theta_m(IJK,M))
708     
709     ! Bulk viscosity
710     ! The following formulation is a guess for Mu_b be by taking 5/3 of the
711     ! collisional viscosity contribution. In this case col. visc. is the
712     ! eps^2 contribution to Mu_s_v(IJK). This might be changed later if
713     ! communications with Ahmadi reveals a different appoach
714                 Mu_b_v(IJK) = 5.d0/3.d0* Tmp_Ahmadi_Const                  &
715                     *0.1045d0*(12.1824d0*G_0(IJK,M,M)*EP_s(IJK,M)*EP_s(IJK,M)) &
716                     *D_p(IJK,M)*RO_S(IJK,M)* DSQRT(Theta_m(IJK,M))
717     
718     ! Second viscosity as defined in MFIX
719                 LAMBDA_S_V(IJK) = Eta*Mu_b_v(IJK) - (2d0*Mu_s_v(IJK))/3d0
720     
721     ! Defining Ahmadi conductivity from his equation 42 in Cao and Ahmadi 1995 paper
722     ! note the constant 0.0711 is now 0.1306 because K = 3/2 theta_m
723                 Kth_s(IJK,M) = 0.1306D0*RO_S(IJK,M)*D_p(IJK,M)*(ONE+C_e**2)* &
724                     (ONE/G_0(IJK,M,M)+4.8D0*EP_s(IJK,M)+12.1184D0 &
725                     *EP_s(IJK,M)*EP_s(IJK,M)*G_0(IJK,M,M) )  &
726                     *DSQRT(Theta_m(IJK,M))
727     
728     
729              ENDIF   ! Fluid_at
730           ENDDO   ! IJK loop
731           RETURN
732           END SUBROUTINE GT_PDE_AHMADI
733     
734     
735     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
736     !                                                                      C
737     !  Subroutine: GT_PDE_simonin                                          C
738     !  Purpose: Calculate granular stress terms (viscosity, bulk viscosity C
739     !  solids pressure) & granular conductivity using simonin model        C
740     !                                                                      C
741     !  Author: Sofiane Benyahia, Fluent Inc.            Date: 02-01-05     C
742     !                                                                      C
743     !  Literature/Document References:                                     C
744     !  Simonin, O., 1996. Combustion and turbulence in two-phase flows,    C
745     !     Von Karman institute for fluid dynamics, lecture series,         C
746     !     1996-02                                                          C
747     !  Balzer, G., Simonin, O., Boelle, A., and Lavieville, J., 1996,      C
748     !     A unifying modelling approach for the numerical prediction       C
749     !     of dilute and dense gas-solid two phase flow. CFB5, 5th int.     C
750     !     conf. on circulating fluidized beds, Beijing, China.             C
751     !                                                                      C
752     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
753           Subroutine gt_pde_simonin (M)
754     
755     ! Modules
756     !---------------------------------------------------------------------//
757           USE constant, only: C_e, eta
758           USE constant, only: pi
759     
760           USE drag, only: f_gs
761     
762           USE fldvar, only: p_s_v
763           USE fldvar, only: ep_g, ro_g
764           USE fldvar, only: ROP_s, ro_s, d_p, theta_m
765           USE fldvar, only: ep_s
766           USE fldvar, only: k_turb_g, e_turb_g
767     
768           USE kintheory, only: kt_rvel, kt_dga, kt_cos_theta
769           USE param1, only: zero, one, large_number, small_number
770     
771           USE physprop, only: smax
772           USE physprop, only: kth_s
773     
774           USE rdf, only: g_0
775     
776           USE toleranc, only: dil_ep_s
777     
778           USE turb, only: tau_1, tau_12, k_12
779     
780           USE visc_s, only: lambda_s_v, mu_s_v, mu_b_v
781     
782           USE functions, only: fluid_at
783           USE compar, only: ijkstart3, ijkend3
784           IMPLICIT NONE
785     
786     ! Dummy arguments
787     !---------------------------------------------------------------------//
788     ! solids phase index
789           INTEGER, INTENT(IN) :: M
790     
791     ! Local parameters
792     !---------------------------------------------------------------------//
793     ! constant in simonin model
794           DOUBLE PRECISION, PARAMETER :: C_mu = 9.0D-02
795     
796     ! Local variables
797     !---------------------------------------------------------------------//
798     ! cell indices
799           INTEGER :: IJK
800     ! solids phase index
801           INTEGER :: MM, L
802     ! defining parameters for Simonin model
803           DOUBLE PRECISION :: Tau_12_st, Tau_2_c, Tau_2, Zeta_r, C_Beta
804           DOUBLE PRECISION :: Sigma_c, Zeta_c, Omega_c, Zeta_c_2, X_21, Nu_t
805           DOUBLE PRECISION :: MU_2_T_Kin, Mu_2_Col, Kappa_kin, Kappa_Col
806           DOUBLE PRECISION :: cos_theta
807     ! sum of ep_s * g_0
808           DOUBLE PRECISION :: SUM_EpsGo
809     ! single particle drag coefficient/ep_s
810           DOUBLE PRECISION :: dga_sl
811     ! relative velocity between gas and solids
812           DOUBLE PRECISION :: rvel_l
813     !---------------------------------------------------------------------//
814     
815           DO IJK = ijkstart3, ijkend3
816              IF ( FLUID_AT(IJK) ) THEN
817     
818     ! Define time scales and constants
819     
820     ! time scale of turbulent eddies
821                 Tau_1(ijk) = 3.d0/2.d0*C_MU*K_Turb_G(IJK)/&
822                    (E_Turb_G(IJK)+small_number)
823     
824     ! NOTE: Some calculations are based explicitly on solids phase 1!
825     ! Parameters based on L=1: tau_12_st, zeta_r, cos_theta, c_beta,
826     ! tau_12, tau_2, nu_t, k_12...
827                 L = 1
828     
829     ! Particle relaxation time. For very dilute flows avoid singularity
830     ! by redefining the drag as single particle drag
831                 IF(Ep_s(IJK,M) > DIL_EP_S .AND. &
832                    F_GS(IJK,L) > small_number) THEN
833                    Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,L)
834                 ELSE
835                    DgA_sL = KT_DGA(ijk, l)
836                    Tau_12_st = RO_S(IJK,M)/DgA_sL
837                 ENDIF
838     
839     ! This is Zeta_r**2 as defined by Simonin
840                 rvel_l = KT_RVEL(ijk, l)
841                 Zeta_r = 3.0d0 * RVEL_L**2 / &
842                    (2.0d0*K_Turb_G(IJK)+small_number)
843     
844                 cos_theta = KT_COS_Theta(IJK, L)
845                 C_Beta = 1.8d0 - 1.35d0*Cos_Theta**2
846     
847     ! Lagrangian Integral time scale: Tau_12
848     ! the time-scale of the fluid turbulent motion viewed by the
849     ! particle (crossing trajectory effect):
850     ! no solids tau_12 = tau_1
851                 Tau_12(ijk) = Tau_1(ijk)/sqrt(ONE+C_Beta*Zeta_r)
852     
853     ! Defining the inter-particle collision time
854                 IF(Ep_s(IJK,M) > DIL_EP_S) THEN
855                    Tau_2_c = D_p(IJK,M)/(6.d0*Ep_s(IJK,M)*G_0(IJK,M,M) &
856                    *DSQRT(16.d0*(Theta_m(ijk,m)+Small_number)/PI))
857                 ELSE             ! assign it a large number
858                    Tau_2_c = LARGE_NUMBER
859                 ENDIF
860     
861     ! Define some constants
862                 Sigma_c = (ONE+ C_e)*(3.d0-C_e)/5.d0
863     ! Zeta_c: const. to be used in the K_2 Diffusion coefficient.
864                 Zeta_c  = (ONE+ C_e)*(49.d0-33.d0*C_e)/100.d0
865                 Omega_c = 3.d0*(ONE+ C_e)**2 *(2.d0*C_e-ONE)/5.d0
866                 Zeta_c_2= 2./5.*(ONE+ C_e)*(3.d0*C_e-ONE)
867     
868     ! mixed time scale in the generalized Simonin theory (switch between dilute
869     ! and kinetic theory formulation of the stresses)
870                 Tau_2 = ONE/(2./Tau_12_st+Sigma_c/Tau_2_c)
871     ! The ratio of these two time scales.
872                 Nu_t =  Tau_12(ijk)/Tau_12_st
873     
874     
875     ! The ratio of densities
876                 X_21 = Ep_s(IJK,M)*RO_S(IJK,M)/(EP_g(IJK)*RO_g(IJK))
877     
878     ! Definition of an "algebraic" form of of Simonin K_12 PDE. This is obtained
879     ! by equating the dissipation term to the exchange terms in the PDE and
880     ! neglecting all other terms, i.e. production, convection and diffusion.
881     ! This works because Tau_12 is very small for heavy particles
882                 K_12(ijk) = Nu_t / (ONE+Nu_t*(ONE+X_21)) * &
883                     (2.d+0 *K_Turb_G(IJK) + 3.d+0 *X_21*theta_m(ijk,m))
884     
885     ! Realizability Criteria
886                 IF(K_12(ijk) > DSQRT(6.0D0*K_Turb_G(IJK)*theta_m(ijk,m))) THEN
887                    K_12(ijk) = DSQRT(6.0D0*K_Turb_G(IJK)*theta_m(ijk,m))
888                 ENDIF
889     
890     ! The following is purely an ad-hoc modification so that the underlying
891     ! monodisperse theory can be used for polydisperse systems in a
892     ! consistent manner.  That is, solids pressure, viscosity and
893     ! conductivity must be additive. The non-linear terms (eps^2) are
894     ! corrected so the stresses of two or more identical solids phases are
895     ! equal to those of a equivalent single solids phase. sof June 15 2005.
896                 SUM_EpsGo = ZERO
897                 DO MM = 1, SMAX
898                    SUM_EpsGo =  SUM_EpsGo+EP_s(IJK,MM)*G_0(IJK,M,MM)
899                 ENDDO
900     
901     
902     ! Solids pressure
903     ! Note this formulation is the same as standard granular pressure
904                 P_s_v(IJK) = ROP_s(IJK,M) * &
905                       (ONE + 4.d0*Eta*SUM_EpsGo)*Theta_m(IJK,M)
906     
907     ! Solids viscosity: shear and bulk
908     ! Turbulent Kinetic (MU_2_T_Kin) and collisional (Mu_2_Col) viscosities
909                 MU_2_T_Kin = (2.0d0/3.0d0*K_12(ijk)*Nu_t + Theta_m(IJK,M) * &
910                     (ONE+ zeta_c_2*EP_s(IJK,M)*G_0(IJK,M,M)))*Tau_2
911                 Mu_2_Col = 8.d0/5.d0*EP_s(IJK,M)*G_0(IJK,M,M)*Eta* (MU_2_T_Kin+ &
912                     D_p(IJK,M)*DSQRT(Theta_m(IJK,M)/PI))
913                 Mu_b_v(IJK) = 5.d0/3.d0*EP_s(IJK,M)*RO_S(IJK,M)*Mu_2_Col
914                 Mu_s_v(IJK) = EP_s(IJK,M)*RO_S(IJK,M)*(MU_2_T_Kin + Mu_2_Col)
915     
916     
917     ! Second viscosity as defined in MFIX
918                 LAMBDA_S_V(IJK) = Eta*Mu_b_v(IJK) - (2d0*Mu_s_v(IJK))/3d0
919     
920     ! Solids conductivity
921     ! Defining Turbulent Kinetic diffusivity: Kappa
922                 Kappa_kin = (9.d0/10.d0*K_12(ijk)*Nu_t + 3.0D0/2.0D0 * &
923                     Theta_m(IJK,M)*(ONE+ Omega_c*EP_s(IJK,M)*G_0(IJK,M,M)))/&
924                     (9.d0/(5.d0*Tau_12_st) + zeta_c/Tau_2_c)
925                 Kappa_Col = 18.d0/5.d0*EP_s(IJK,M)*G_0(IJK,M,M)*Eta* &
926                     (Kappa_kin+ 5.d0/9.d0*D_p(IJK,M)*DSQRT(Theta_m(IJK,M)/PI))
927                 Kth_s(IJK,M) =  EP_s(IJK,M)*RO_S(IJK,M)*(Kappa_kin + Kappa_Col)
928     
929     
930              ENDIF   ! Fluid_at
931           ENDDO   ! IJK loop
932           RETURN
933           END SUBROUTINE GT_PDE_SIMONIN
934     
935     
936     
937     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
938     !                                                                      C
939     !  Subroutine: GT_PDE_GD                                               C
940     !  Purpose: Implement kinetic theory of Garzo and Dufty (1999) for     C
941     !  calculation of granular stress terms and granular conductivity      C
942     !                                                                      C
943     !  Author: Janine E. Galvin                                            C
944     !                                                                      C
945     !  Literature/Document References:                                     C
946     !  Garzo, V., and Dufty, J., Homogeneous cooling state for a           C
947     !     granular mixture, Physical Review E, 1999, Vol 60 (5), 5706-     C
948     !     5713                                                             C
949     !                                                                      C
950     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
951           Subroutine gt_pde_gd (M)
952     
953     ! Modules
954     !---------------------------------------------------------------------//
955           USE constant, only: C_e, pi, switch
956           USE drag, only: f_gs
957           USE fldvar, only: p_s_v
958           USE fldvar, only: rop_s, ro_s, d_p, theta_m
959           USE fldvar, only: ep_s
960           USE fldvar, only: ro_g
961           USE kintheory, only: kt_dga
962           USE param1, only: zero, small_number
963           USE physprop, only: kth_s, kphi_s
964           USE rdf, only: g_0, dg_0dnu
965           USE toleranc, only: dil_ep_s
966           USE visc_s, only: lambda_s_v, mu_s_v, mu_b_v
967           USE compar, only: ijkstart3, ijkend3
968           USE functions, only: fluid_at
969           IMPLICIT NONE
970     
971     ! Dummy arguments
972     !---------------------------------------------------------------------//
973     ! solids phase index
974           INTEGER, INTENT(IN) :: M
975     
976     ! Local variables
977     !---------------------------------------------------------------------//
978     ! cell Indices
979           INTEGER :: IJK
980     ! Use to compute MU_s(IJK,M) & Kth_S(IJK,M)
981           DOUBLE PRECISION :: Mu_star, Kth_star
982     !
983           DOUBLE PRECISION :: D_PM, M_PM, NU_PM, EP_SM, RO_SM, ROP_SM
984     !
985           DOUBLE PRECISION :: chi, dChiOdphi
986     !
987           DOUBLE PRECISION :: c_star, zeta0_star, nu_eta_star, &
988                               gamma_star, eta_k_star, eta_star, eta0, &
989                               kappa0, nu_kappa_star, kappa_k_star, &
990                               qmu_k_star, qmu_star, kappa_star, press_star
991     ! single particle drag coefficient
992           DOUBLE PRECISION :: dga_sm
993     !---------------------------------------------------------------------//
994     
995           DO IJK = ijkstart3, ijkend3
996               IF ( FLUID_AT(IJK) ) THEN
997     
998     ! local aliases
999                  D_PM = D_P(IJK,M)
1000                  RO_SM = RO_S(IJK,M)
1001                  ROP_SM = ROP_S(IJK,M)
1002                  EP_SM = EP_S(IJK,M)
1003                  M_PM = (PI/6.d0)*D_PM**3 * RO_SM
1004                  NU_PM = ROP_SM/M_PM
1005                  Chi = G_0(IJK,M,M)
1006                  dChiOdphi = DG_0DNU(EP_SM)
1007                  IF(EP_SM <= DIL_EP_S) &
1008                     dga_sm = KT_DGA(IJK, M)
1009     
1010     
1011     ! Pressure/Viscosity/Bulk Viscosity
1012     ! Note: k_boltz = M_PM
1013     !-----------------------------------
1014     ! Find pressure in the Mth solids phase
1015                  press_star = 1.d0 + 2.d0*(1.d0+C_E)*EP_SM*Chi
1016     
1017     ! n*k_boltz = n*m = ep_s*ro_s
1018                  P_s_v(IJK) = ROP_sM*Theta_m(IJK,M)*press_star
1019     
1020     ! find bulk and shear viscosity
1021                  c_star = 32.0d0*(1.0d0 - C_E)*(1.d0 - 2.0d0*C_E*C_E) &
1022                     / (81.d0 - 17.d0*C_E + 30.d0*C_E*C_E*(1.0d0-C_E))
1023     
1024                  zeta0_star = (5.d0/12.d0)*Chi*(1.d0 - C_E*C_E) &
1025                     * (1.d0 + (3.d0/32.d0)*c_star)
1026     
1027                  nu_eta_star = Chi*(1.d0 - 0.25d0*(1.d0-C_E)*(1.d0-C_E)) &
1028                     * (1.d0-(c_star/64.d0))
1029     
1030                  gamma_star = (4.d0/5.d0)*(32.d0/PI)*EP_SM*EP_SM &
1031                     * Chi*(1.d0+C_E) * (1.d0 - (c_star/32.d0))
1032     
1033                  eta_k_star = (1.d0 - (2.d0/5.d0)*(1.d0+C_E)*(1.d0-3.d0*C_E) &
1034                     * EP_SM*Chi ) / (nu_eta_star - 0.5d0*zeta0_star)
1035     
1036                  eta_star = eta_k_star*(1.d0 + (4.d0/5.d0)*EP_SM*Chi &
1037                     * (1.d0+C_E) ) + (3.d0/5.d0)*gamma_star
1038     
1039                  eta0 = 5.0d0*M_PM*DSQRT(Theta_m(IJK,M)/PI) / (16.d0*D_PM*D_PM)
1040     
1041                  IF(SWITCH == ZERO .OR. RO_G(IJK) == ZERO) THEN
1042                     Mu_star = eta0
1043                  ELSEIF(Theta_m(IJK,M) .LT. SMALL_NUMBER)THEN
1044                     Mu_star = ZERO
1045                  ELSEIF(EP_SM <= DIL_EP_S) THEN
1046                     Mu_star = RO_SM*EP_SM*Chi*Theta_m(IJK,M)*eta0 / &
1047                        ( RO_S(IJK,M)*EP_SM*Chi*Theta_m(IJK,M) + &
1048                        2.d0*DgA_sM*eta0/RO_S(IJK,M) )
1049                  ELSE
1050                     Mu_star = RO_SM*EP_SM*Chi*Theta_m(IJK,M)*eta0 / &
1051                        ( RO_SM*EP_SM*Chi*Theta_m(IJK,M) + &
1052                        (2.d0*F_gs(IJK,M)*eta0/(RO_SM*EP_SM)) )
1053                  ENDIF
1054     
1055     ! Shear and bulk viscosity
1056                  Mu_s_v(IJK) = Mu_star * eta_star
1057                  Mu_b_v(IJK) = Mu_star * gamma_star
1058     
1059     ! Second viscosity
1060                  LAMBDA_S_V(IJK) = Mu_b_v(IJK) - (2.d0/3.d0)*Mu_s_v(IJK)
1061     
1062     
1063     ! Granular Conductivity/Dufour Coefficient
1064     !-----------------------------------
1065                  kappa0 = (15.d0/4.d0)*eta0
1066     
1067                  nu_kappa_star = (Chi/3.d0)*(1.d0+C_E) * ( 1.d0 + &
1068                     (33.d0/16.d0)*(1.d0-C_E) + ((19.d0-3.d0*C_E)/1024.d0)*&
1069                     c_star)
1070     !             nu_mu_star = nu_kappa_star
1071     
1072                  kappa_k_star = (2.d0/3.d0)*(1.d0 +0.5d0*(1.d0+press_star)*&
1073                     c_star + (3.d0/5.d0)*EP_SM*Chi*(1.d0+C_E)*(1.d0+C_E) * &
1074                     (2.d0*C_E - 1.d0 + ( 0.5d0*(1.d0+C_E) - 5.d0/&
1075                     (3*(1.d0+C_E))) * c_star ) ) / (nu_kappa_star - &
1076                     2.d0*zeta0_star)
1077     
1078                  kappa_star = kappa_k_star * (1.d0 + (6.d0/5.d0)*EP_SM* &
1079                     Chi*(1.d0+C_E) ) + (256.d0/25.d0)*(EP_SM* &
1080                     EP_SM/PI)*Chi*(1.d0+C_E)*(1.d0+(7.d0/32.d0)* &
1081                     c_star)
1082     
1083                  IF(SWITCH == ZERO .OR. RO_G(IJK) == ZERO) THEN
1084                     Kth_star= kappa0
1085                  ELSEIF(Theta_m(IJK,M) .LT. SMALL_NUMBER)THEN
1086                     Kth_star = ZERO
1087                  ELSEIF(EP_SM <= DIL_EP_S) THEN
1088                     Kth_star = RO_SM*EP_SM*Chi*Theta_m(IJK,M)*kappa0/ &
1089                        (RO_SM*EP_SM*Chi*Theta_m(IJK,M) + 1.2d0*DgA_sM* &
1090                        kappa0/RO_SM)
1091                  ELSE
1092                     Kth_star = RO_SM*EP_SM*Chi*Theta_m(IJK,M)*kappa0/ &
1093                        (RO_SM*EP_SM*Chi*Theta_m(IJK,M)+ &
1094                        (1.2d0*F_gs(IJK,M)*kappa0/(RO_SM*EP_SM)) )
1095                  ENDIF
1096     
1097     ! Granular conductivity
1098                  Kth_s(IJK,M) = Kth_star * kappa_star
1099     
1100     ! transport coefficient of the Mth solids phase associated
1101     ! with gradient in volume fraction in heat flux
1102                  qmu_k_star = 2.d0*( (1.d0+EP_SM*dChiOdphi)* &
1103                     zeta0_star*kappa_k_star + ( (press_star/3.d0) + &
1104                     (2.d0/3.d0)* EP_SM*(1.d0+C_E)*(Chi+EP_SM* dChiOdphi) )*&
1105                     c_star - (4.d0/5.d0)*EP_SM*Chi* (1.d0+(EP_SM/2.d0)*&
1106                     dChiOdphi)* (1.d0+C_E) * ( C_E*(1.d0-C_E)+0.25d0*&
1107                     ((4.d0/3.d0)+C_E* (1.d0-C_E))*c_star ) ) / &
1108                     (2.d0*nu_kappa_star-3.d0*zeta0_star)
1109     
1110                  qmu_star = qmu_k_star*(1.d0+(6.d0/5.d0)*EP_SM*Chi*&
1111                     (1.d0+C_E) )
1112     
1113                  IF (EP_SM .LT. SMALL_NUMBER) THEN
1114                     Kphi_s(IJK,M) = ZERO
1115                  ELSE
1116                     Kphi_s(IJK,M) = (Theta_m(IJK,M)*Kth_star/NU_PM)*qmu_star
1117                  ENDIF
1118     
1119               ENDIF   ! Fluid_at
1120           ENDDO   ! IJK loop
1121           RETURN
1122           END SUBROUTINE GT_PDE_GD
1123     
1124     
1125     
1126     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1127     !                                                                      C
1128     !  Subroutine: GT_PDE_GTSH                                             C
1129     !  Purpose: Implement kinetic theory of Garzo, Tenneti, Subramaniam    C
1130     !  Hrenya (2012) for calculation of granular stress terms and          C
1131     !                                                                      C
1132     !  Author: Sofiane Benyahia                                            C
1133     !                                                                      C
1134     !  Literature/Document References:                                     C
1135     !  Garzo, V., Tenneti, S., Subramaniam, S., and Hrenya, C. M.,         C
1136     !      "Enskog kinetic theory for monodisperse gas-solid flows", JFM,  C
1137     !      Vol. 712, 2012, pp. 129-168                                     C
1138     !                                                                      C
1139     !  Comments:                                                           C
1140     !  And also based on C.M. Hrenya hand-notes dated Sep 2013             C
1141     !                                                                      C
1142     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1143           Subroutine gt_pde_gtsh (M)
1144     
1145     ! Modules
1146     !---------------------------------------------------------------------//
1147           USE constant, only: C_e, pi
1148     
1149           USE fldvar, only: p_s_v
1150           USE fldvar, only: ro_g
1151           USE fldvar, only: rop_s, ro_s, d_p, theta_m
1152           USE fldvar, only: ep_s
1153     
1154           USE kintheory, only: EDT_s_ip, xsi_gtsh, a2_gtsh
1155           USE kintheory, only: kt_rvel
1156           USE kintheory, only: epm, G_gtsh, S_star, K_phi, R_d
1157     
1158           USE param1, only: zero, one, small_number
1159     
1160           USE physprop, only: mu_g
1161           USE physprop, only: kth_s, kphi_s
1162     
1163           USE rdf, only: g_0, dg_0dnu
1164     
1165           USE visc_s, only: lambda_s_v, mu_s_v, mu_b_v
1166     
1167           USE compar, only: ijkstart3, ijkend3
1168           USE functions, only: fluid_at
1169           IMPLICIT NONE
1170     
1171     ! Dummy arguments
1172     !---------------------------------------------------------------------//
1173     ! solids phase index
1174           INTEGER, INTENT(IN) :: M
1175     
1176     ! Local variables
1177     !---------------------------------------------------------------------//
1178     ! cell Indices
1179           INTEGER :: IJK
1180     !
1181           DOUBLE PRECISION :: D_PM, M_PM, NU_PM, EP_SM
1182     !
1183           DOUBLE PRECISION :: chi, dChiOdphi, eta0
1184     !
1185           DOUBLE PRECISION :: nu0, nuN, etaK
1186           DOUBLE PRECISION :: dZeta_dT, dGama_dT, NuK, Kth0, KthK
1187           DOUBLE PRECISION :: Rdissdphi, Kphidphi, Re_T, dGamadn, dRdphi
1188           DOUBLE PRECISION :: denom
1189           DOUBLE PRECISION :: dSdphi, R_dphi, Tau_st, dPsidn, MuK
1190     ! relative velocity
1191           DOUBLE PRECISION :: RVEL
1192     !---------------------------------------------------------------------//
1193     
1194           DO IJK = ijkstart3, ijkend3
1195               IF ( FLUID_AT(IJK) ) THEN
1196     
1197     ! local aliases
1198                  D_PM = D_P(IJK,M)
1199                  EP_SM = EP_S(IJK,M)
1200                  M_PM = (PI/6.d0)*D_PM**3 * RO_S(IJK,M)
1201                  NU_PM = ROP_S(IJK,M)/M_PM
1202                  Chi = G_0(IJK,M,M)
1203                  dChiOdphi = DG_0DNU(EP_SM)
1204                  RVEL = kt_rvel(ijk, m)
1205     
1206     
1207     ! note that T = (m_pm*theta_m)
1208     ! Pressure/Viscosity/Bulk Viscosity
1209     !-----------------------------------
1210     ! solids pressure, eq (6.14) of GTSH theory
1211                  P_s_v(IJK) = ROP_s(IJK,M)*Theta_m(IJK,M)* &
1212                     (one+2d0*(one+C_E)*Chi*EP_SM)
1213     
1214     ! evaluating shear viscosity, eq (7.3) of GTSH theory
1215     ! starting with nu_0 equ (7.6-7.8)
1216                  eta0 = 0.3125d0/(dsqrt(pi)*D_PM**2)*M_pm*&
1217                     dsqrt(theta_m(ijk,m))
1218     
1219                  nu0 = (96.d0/5.d0)*(EP_SM/D_PM)*DSQRT(Theta_m(IJK,M)/PI)
1220     !             nu0 = NU_PM*M_pm*theta_m(ijk,m)/eta0
1221     
1222                  nuN = 0.25d0*nu0*Chi*(3d0-C_E)*(one+C_E) * &
1223                         (one+0.4375d0*A2_gtsh(ijk))
1224     ! defining kinetic part of shear viscosity nuK  equ (7.7)
1225                  etaK = rop_s(ijk,m)*theta_m(ijk,m) / (nuN-0.5d0*( &
1226                     EDT_s_ip(ijk,M,M)-xsi_gtsh(ijk)/theta_m(ijk,m) - &
1227                     2d0*G_gtsh(EP_SM, Chi, IJK, M)/M_PM)) * (one -0.4d0 * &
1228                     (one+C_E)*(one-3d0*C_E)*EP_SM*Chi)
1229     
1230     ! bulk viscosity lambda eq. (7.5)
1231                  Mu_b_v(IJK) = 25.6d0/pi * EP_SM**2 * Chi *(one+C_E) * &
1232                     (one - A2_gtsh(ijk)/16d0)*eta0
1233     
1234     ! Finally shear viscosity, eq (7.9) of GTSH theory
1235                  Mu_s_v(IJK) = etaK*(one+0.8d0*EP_SM*Chi*(one+C_E)) + &
1236                     0.6d0*Mu_b_v(IJK)
1237     
1238     ! Second viscosity as defined in MFIX
1239                  LAMBDA_S_V(IJK) = Mu_b_v(IJK) - (2.d0/3.d0)*Mu_s_v(IJK)
1240     
1241     
1242     ! Conductivity/Dufour coefficient
1243     !-----------------------------------
1244     ! Calculate conductivity Kth_s(IJK,M), eq. 7.12 GTSH theory.
1245     ! Start with calculating dZeta/dT and dGama/dT
1246     ! note that 1/Tau**2 = (3d0*pi*mu_g(ijk)*D_PM/M_p)**2 defined
1247     ! under eq. 8.2 GTSH
1248                  dZeta_dT = -0.5d0*xsi_gtsh(ijk)/(M_pm*theta_m(ijk,m))
1249     
1250                  dGama_dT = 3d0*pi*D_PM**2*RO_g(ijk)*K_phi(EP_SM)/ &
1251                     (2d0*M_pm*dsqrt(theta_m(ijk,m)))
1252     
1253     ! evaluating eq (7.16) in GTSH
1254                  NuK = nu0*(one+C_E)/3d0*Chi*( one+2.0625d0*(one-C_E)+ &
1255                     ((947d0-579*C_E)/256d0*A2_gtsh(ijk)) )
1256     
1257     ! evaluating eq. (7.13)
1258                  Kth0 = 3.75d0*eta0/M_pm
1259     
1260     ! evaluating kinetic conductivity Kk eq. (7.14)
1261     ! note that 1/2m/T Psi and m dZeta_dT cancel out.
1262                  KthK = zero
1263                  IF(EP_SM > SMALL_NUMBER) KthK = 2d0/3d0*Kth0*nu0 / (NuK - &
1264                     2d0*EDT_s_ip(ijk,M,M) - 2d0*theta_m(ijk,m)*dGama_dT) * &
1265                     (one+2d0*A2_gtsh(ijk)+0.6d0*EP_SM*Chi* &
1266                     (one+C_E)**2*(2*C_E-one+A2_gtsh(ijk)*(one+C_E)))
1267     
1268     ! the conductivity Kth from eq (7.17) in GTSH theory:
1269                  Kth_s(IJK,M) = KthK*(one+1.2d0*EP_SM*Chi*(one+C_E)) + &
1270                     (10.24d0/pi* EP_SM**2*Chi*(one+C_E)*(one+0.4375d0* &
1271                     A2_gtsh(ijk))*Kth0)
1272     
1273     ! Finaly notice that conductivity K in eq (7.10) must be
1274     ! multiplied by m because of grad(T)
1275                  Kth_s(IJK,M) = M_pm * Kth_s(IJK,M)
1276     
1277     ! Calculate the Dufour coefficient Kphi_s(IJK,M) in equation (7.18) of
1278     ! GTSH theory.
1279     ! First, calculate terms in 2 n/m x Gama_n, dRdiss/dphi and dK_phi/dphi.
1280     ! Notice that 2 n/m Gama_n = 2 phi/m Gama_phi, so multiply the
1281     ! deriviatives of Rdiss and K_phi by phi to avoid possible division by
1282     ! phi.
1283                  Rdissdphi = ZERO
1284                  IF(EP_SM > SMALL_NUMBER) Rdissdphi = &
1285                     1.5d0*dsqrt(EP_SM/2d0)+135d0/64d0*EP_SM*(dlog(EP_SM)+one) +&
1286                     11.26d0*EP_SM*(one-10.2*EP_SM+49.71d0*EP_SM**2-87.08d0* &
1287                     EP_SM**3) - EP_SM*dlog(epM)*(Chi+EP_SM*dChiOdphi)
1288     ! corrections due to W. Fullmer
1289                  Kphidphi = EP_SM*(0.212d0*0.142d0/(EP_SM**0.788d0*&
1290                     (one-EP_SM)**4.454d0) + 4.454d0*K_phi(EP_SM)/(one-EP_SM))
1291                  Kphidphi = zero  ! this is compatible with K_phi = zero
1292     
1293                  Re_T = ro_g(ijk)*D_PM*dsqrt(theta_m(ijk,m)) / mu_g(ijk)
1294     
1295     ! The term phi x Gama_phi becomes
1296                  dGamadn = 3d0*pi*D_pm*Mu_g(ijk)*(Rdissdphi+Re_T*Kphidphi)
1297     
1298     ! Second, calculate terms in 2 rho x Psi_n, which is same as ro_s x
1299     ! EP_SM x Psi_n
1300     ! Take EP_SM inside the derivative dS_star/dphi to avoid singularities
1301     ! as EP_SM -> 0
1302     
1303     ! calculating the term phi*dRd/dphi
1304                  dRdphi = zero
1305                  IF((EP_SM > SMALL_NUMBER) .AND. (EP_SM <= 0.4d0)) THEN
1306                     denom = one+0.681d0*EP_SM-8.48d0*EP_SM**2+8.16d0*EP_SM**3
1307                     dRdphi = (1.5d0*dsqrt(EP_SM/2d0)+135d0/64d0*EP_SM*&
1308                        (dlog(EP_SM)+one)+ 17.14d0*EP_SM)/denom - EP_SM*&
1309                        (one+3d0*dsqrt(EP_SM/2d0) + 135d0/64d0*EP_SM*&
1310                        dlog(EP_SM)+17.14*EP_SM)/denom**2 * &
1311                        (0.681d0-16.96d0*EP_SM+24.48d0*EP_SM**2)
1312                  ELSEIF(EP_SM > 0.4d0) THEN
1313                     dRdphi = 10d0*(one+2d0*EP_SM)/(one-EP_SM)**4
1314                  ENDIF
1315     
1316     ! calculating the term phi*dS_star/dphi
1317                  dSdphi = zero
1318                  IF(EP_SM >= 0.1d0) THEN
1319                     R_dphi = R_d(EP_SM)
1320                     denom = one+3.5d0*dsqrt(EP_SM)+5.9d0*EP_SM
1321                     dSdphi = 2d0*R_dphi*dRdphi/(Chi*denom) - &
1322                        EP_SM*R_dphi**2 * (dChiOdphi/(Chi**2*denom) + &
1323                        (1.75d0/dsqrt(EP_SM)+5.9d0)/(Chi*denom**2))
1324                  ENDIF
1325     
1326     ! defining the relaxation time Tau_st
1327                  Tau_st = M_pm/(3d0*pi*mu_g(ijk)*D_pm)
1328     
1329     ! The term phi x Psi_n becomes
1330                  dPsidn = dsqrt(pi)*D_pm**4*RVEL**2 / &
1331                     (36d0*Tau_st**2*dsqrt(theta_m(ijk,m))) * dSdphi
1332     
1333     ! Now compute the kinetic contribution to Dufour coef. Muk eq (7.20) GSTH
1334                  Muk = ZERO  ! This is assumed to avoid /0 for EP_SM = 0
1335                  IF(EP_SM > SMALL_NUMBER) Muk = Kth0*Nu0*M_pm*&
1336                     theta_m(ijk,m)/NU_PM / (NuK-1.5d0*(EDT_s_ip(ijk,M,M)-&
1337                     xsi_gtsh(ijk)/theta_m(ijk,m))) * ( KthK/(Kth0*Nu0)*&
1338                     (2d0/M_pm*dGamadn-ro_s(IJK,m)/(M_pm* theta_m(ijk,m))*&
1339                     dPsidn + EDT_s_ip(ijk,M,M)*(one+EP_SM/Chi*dChiOdphi)) +&
1340                     2d0/3d0*A2_gtsh(ijk) + 0.8d0*EP_SM*Chi* (one+C_E)*&
1341                     (one+0.5d0*EP_SM/Chi*dChiOdphi)*(C_E*(C_E-one)+ &
1342                     A2_gtsh(ijk)/6d0*(16d0-3d0*C_E+3d0*C_E**2)))
1343     
1344     ! Finaly compute the Dufour coefficient Mu (Kphi_s(IJK,M)) from eq (7.22) GTSH
1345                  Kphi_s(IJK,M) = Muk*(one+1.2d0*EP_SM*Chi*(one+C_E))
1346     
1347               ENDIF   ! Fluid_at
1348           ENDDO   ! outer IJK loop
1349     
1350           RETURN
1351           END SUBROUTINE GT_PDE_GTSH
1352     
1353     
1354     
1355     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1356     !                                                                      C
1357     !  Subroutine: GT_PDE_IA                                               C
1358     !  Purpose: Implement kinetic theory of Iddir & Arastoopour (2005) for C
1359     !  calculation of granular stress terms and granular conductivity      C
1360     !                                                                      C
1361     !  Author: Janine E. Galvin, Univeristy of Colorado                    C
1362     !                                                                      C
1363     !  Literature/Document References:                                     C
1364     !  Iddir, Y.H., PhD Modeling of the multiphase mixture of particles    C
1365     !     using the kinetic theory approach, PhD Dissertation in           C
1366     !     Chemical Engineering, Illinois Institute of Technology, 2004     C
1367     !  Iddir, Y.H., H. Arastoopour, and C.M. Hrenya, Analysis of binary    C
1368     !     and ternary granular mixtures behavior using the kinetic         C
1369     !     theory approach. Powder Technology, 2005, p. 117-125.            C
1370     !                                                                      C
1371     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1372           Subroutine gt_pde_ia (M)
1373     
1374     ! Modules
1375     !---------------------------------------------------------------------//
1376           USE constant, only: switch, switch_ia
1377           USE constant, only: C_e
1378           USE constant, only: pi
1379     
1380           USE drag, only: f_gs
1381           USE fldvar, only: RO_g
1382           USE fldvar, only: ROP_s, RO_s, D_p, theta_m
1383           USE fldvar, only: p_s_v
1384           USE fldvar, only: ep_s
1385     
1386           USE kintheory, only: mu_sm_ip, mu_sl_ip
1387           USE kintheory, only: xi_sm_ip, xi_sl_ip
1388           USE kintheory, only: kth_sl_ip, knu_sm_ip, knu_sl_ip, kvel_s_ip
1389           USE kintheory, only: kt_dga
1390     
1391           USE param1, only: one, zero, small_number
1392     
1393           USE physprop, only: smax
1394           USE physprop, only: Kth_s
1395     
1396           USE rdf, only: g_0
1397     
1398           USE toleranc, only: dil_ep_s
1399           USE ur_facs, only: UR_Kth_sml
1400           USE visc_s, only: mu_s_v, lambda_s_v
1401     
1402           USE compar, only: ijkstart3, ijkend3
1403           USE functions, only: fluid_at
1404           IMPLICIT NONE
1405     
1406     ! Dummy arguments
1407     !---------------------------------------------------------------------//
1408     ! solids phase index
1409           INTEGER, INTENT(IN) :: M
1410     
1411     ! Local variables
1412     !---------------------------------------------------------------------//
1413     ! cell indices
1414           INTEGER :: IJK
1415     ! solids phase index
1416           INTEGER :: L
1417     ! use to compute MU_s(IJK,M) & Kth_S(IJK,M)
1418           DOUBLE PRECISION :: Mu_star, Mu_s_dil, Kth_star, K_s_dil
1419     ! variables for Iddir equipartition model
1420           DOUBLE PRECISION :: P_s_sum, P_s_MM, P_s_LM
1421           DOUBLE PRECISION :: MU_common_term, K_common_term
1422           DOUBLE PRECISION :: Mu_sM_sum, MU_s_MM, MU_s_LM, MU_sM_LM, MU_sL_LM
1423           DOUBLE PRECISION :: XI_sM_sum, XI_s_v
1424           DOUBLE PRECISION :: M_PM, M_PL, MPSUM, NU_PL, NU_PM, D_PM, D_PL, DPSUMo2
1425           DOUBLE PRECISION :: Ap_lm, Dp_lm, R0p_lm, R1p_lm, R8p_lm, R9p_lm, Bp_lm,&
1426                               R5p_lm, R6p_lm, R7p_lm
1427           DOUBLE PRECISION :: K_s_sum, K_s_MM, K_s_LM
1428     ! Sum of ep_s * g_0
1429           DOUBLE PRECISION :: SUM_EpsGo
1430     ! Current value of Kth_sl_ip (i.e., without underrelaxation)
1431           DOUBLE PRECISION :: Kth_sL_iptmp
1432     ! drag coefficient/ep_s
1433           DOUBLE PRECISION :: dga_sm
1434     !---------------------------------------------------------------------//
1435     
1436           DO IJK = ijkstart3, ijkend3
1437              IF ( FLUID_AT(IJK) ) THEN
1438     
1439     ! Added for consistency of IA KT (see constant_mod for details)
1440                 IF(SWITCH_IA) THEN
1441                    SUM_EpsGo = ZERO
1442                    DO L = 1, SMAX
1443                       SUM_EpsGo =  SUM_EpsGo+EP_s(IJK,L)*G_0(IJK,M,L)
1444                    ENDDO
1445                 ELSE
1446                    SUM_EpsGo =  EP_s(IJK,M)*G_0(IJK,M,M)
1447                 ENDIF
1448     
1449                 P_s_sum = ZERO
1450                 Mu_sM_sum = ZERO
1451                 XI_sM_sum = ZERO
1452                 IF(EP_S(IJK,M) <= DIL_EP_s) &
1453                    dga_sm = KT_DGA(ijk, M)
1454     
1455                 D_PM = D_P(IJK,M)
1456                 M_PM = (PI/6.d0)*D_PM**3 * RO_S(IJK,M)
1457                 NU_PM = ROP_S(IJK,M)/M_PM
1458     
1459                 P_s_MM = NU_PM*Theta_m(IJK,M)
1460     
1461                 MU_s_dil = (5.d0/96.d0)*D_PM* RO_S(IJK,M)*&
1462                     DSQRT(PI*Theta_m(IJK,M)/M_PM)
1463     
1464                 IF(SWITCH == ZERO .OR. RO_G(IJK) == ZERO) THEN
1465                    Mu_star = MU_s_dil
1466                 ELSEIF(Theta_m(IJK,M)/M_PM < SMALL_NUMBER)THEN
1467                    Mu_star = ZERO
1468                 ELSEIF(EP_S(IJK,M) <= DIL_EP_s) THEN
1469                    Mu_star = MU_s_dil*EP_s(IJK,M)*G_0(IJK,M,M)/ &
1470                       ( SUM_EpsGo + 2.0d0*DgA_sM*MU_s_dil /&
1471                         (RO_S(IJK,M)**2 *(Theta_m(IJK,M)/M_PM)) )
1472                 ELSE
1473                    Mu_star = MU_s_dil*EP_S(IJK,M)*G_0(IJK,M,M)/ &
1474                       ( SUM_EpsGo + 2.0d0*F_gs(IJK,M)*MU_s_dil /&
1475                         (RO_S(IJK,M)**2 *EP_s(IJK,M)*&
1476                          (Theta_m(IJK,M)/M_PM)) )
1477                 ENDIF
1478     
1479                 MU_s_MM = (Mu_star/G_0(IJK,M,M))*&
1480                     (1.d0+(4.d0/5.d0)*(1.d0+C_E)*SUM_EpsGo)**2
1481     
1482                 DO L = 1, SMAX
1483                    D_PL = D_P(IJK,L)
1484                    M_PL = (PI/6.d0)*D_PL**3 * RO_S(IJK,L)
1485                    MPSUM = M_PM + M_PL
1486                    DPSUMo2 = (D_PM+D_PL)/2.d0
1487                    NU_PL = ROP_S(IJK,L)/M_PL
1488     
1489                    IF ( L .eq. M) THEN
1490                       Ap_lm = MPSUM/(2.d0)
1491                       Dp_lm = M_PL*M_PM/(2.d0*MPSUM)
1492                       R0p_lm = ONE/( Ap_lm**1.5 * Dp_lm**2.5 )
1493                       R1p_lm = ONE/( Ap_lm**1.5 * Dp_lm**3 )
1494     
1495                       P_s_LM = PI*(DPSUMo2**3 / 48.d0)*G_0(IJK,M,L)*&
1496                           (M_PM*M_PL/MPSUM)* (M_PM*M_PL)**1.5 *&
1497                           NU_PM*NU_PL*(1.d0+C_E)*R0p_lm*Theta_m(IJK,M)
1498     
1499                       MU_s_LM = DSQRT(PI)*( DPSUMo2**4 / 240d0 )*&
1500                           G_0(IJK,M,L)*(M_PL*M_PM/MPSUM)**2 *&
1501                           (M_PL*M_PM)**1.5 * NU_PM*NU_PL*&
1502                           (1.d0+C_E) * R1p_lm * DSQRT(Theta_m(IJK,M))
1503     
1504     ! This is Mu_i_1 as defined in eq (16) of Galvin document
1505                       MU_sM_ip(IJK,M,L) = (MU_s_MM + MU_s_LM)
1506     
1507     ! This is Mu_ii_2 as defined in eq (17) of Galvin document
1508                       MU_sL_ip(IJK,M,L) = MU_s_LM
1509     
1510     ! solids phase viscosity associated with the trace of
1511     ! solids phase M (eq. 18 from Galvin theory document)
1512                       XI_sM_ip(IJK,M,L) = (5.d0/3.d0)*MU_s_LM
1513     
1514     ! solids phase viscosity associated with the trace
1515     ! of (sum of) all solids phases (eq. 19)
1516                       XI_sL_ip(IJK,M,L) = (5.d0/3.d0)*MU_s_LM
1517     
1518                    ELSE
1519                       Ap_lm = (M_PM*Theta_m(IJK,L)+M_PL*&
1520                           Theta_m(IJK,M))/2.d0
1521                       Bp_lm = (M_PM*M_PL*(Theta_m(IJK,L)-&
1522                           Theta_m(IJK,M) ))/(2.d0*MPSUM)
1523                       Dp_lm = (M_PL*M_PM*(M_PM*Theta_m(IJK,M)+&
1524                           M_PL*Theta_m(IJK,L) ))/&
1525                           (2.d0*MPSUM*MPSUM)
1526                       R0p_lm = (1.d0/(Ap_lm**1.5 * Dp_lm**2.5))+ &
1527                           ((15.d0*Bp_lm*Bp_lm)/(2.d0* Ap_lm**2.5 *&
1528                           Dp_lm**3.5))+&
1529                           ((175.d0*(Bp_lm**4))/(8.d0*Ap_lm**3.5 * &
1530                           Dp_lm**4.5))
1531                       R1p_lm = (1.d0/((Ap_lm**1.5)*(Dp_lm**3)))+ &
1532                           ((9.d0*Bp_lm*Bp_lm)/( Ap_lm**2.5 * Dp_lm**4))+&
1533                           ((30.d0*Bp_lm**4) /( 2.d0*Ap_lm**3.5 * &
1534                           Dp_lm**5))
1535     
1536                       P_s_LM = PI*(DPSUMo2**3 / 48.d0)*G_0(IJK,M,L)*&
1537                           (M_PM*M_PL/MPSUM)* (M_PM*M_PL)**1.5 *&
1538                           NU_PM*NU_PL*(1.d0+C_E)*R0p_lm* &
1539                           (Theta_m(IJK,M)*Theta_m(IJK,L))**2.5
1540     
1541                       MU_common_term = DSQRT(PI)*( DPSUMo2**4 / 240d0 )*&
1542                           G_0(IJK,M,L)*(M_PL*M_PM/MPSUM)**2 *&
1543                           (M_PL*M_PM)**1.5 * NU_PM*NU_PL*&
1544                            (1.d0+C_E) * R1p_lm
1545                       MU_sM_LM = MU_common_term * Theta_m(IJK,M)**2 *&
1546                           Theta_m(IJK,L)**3
1547                       MU_sL_LM = MU_common_term * Theta_m(IJK,L)**2 *&
1548                           Theta_m(IJK,M)**3
1549     
1550     ! solids phase 'viscosity' associated with the divergence
1551     ! of solids phase M. defined in eq (16) of Galvin document
1552                       MU_sM_ip(IJK,M,L) = MU_sM_LM
1553     
1554     ! solids phase 'viscosity' associated with the divergence
1555     ! of all solids phases. defined in eq (17) of Galvin document
1556                       MU_sL_ip(IJK,M,L) = MU_sL_LM
1557     
1558     ! solids phase viscosity associated with the trace of
1559     ! solids phase M
1560                       XI_sM_ip(IJK,M,L) = (5.d0/3.d0)*MU_sM_LM
1561     
1562     ! solids phase viscosity associated with the trace
1563     ! of all solids phases
1564                       XI_sL_ip(IJK,M,L) = (5.d0/3.d0)*MU_sL_LM
1565                    ENDIF
1566     
1567                    P_s_sum = P_s_sum + P_s_LM
1568                    MU_sM_sum = MU_sM_sum + MU_sM_ip(IJK,M,L)
1569                    XI_sM_sum = XI_sM_sum + XI_sM_ip(IJK,M,L)
1570                 ENDDO
1571     
1572     ! Find the term proportional to the identity matrix
1573     ! (pressure in the Mth solids phase)
1574                 P_s_v(IJK) = P_s_sum + P_S_MM
1575     
1576     ! Find the term proportional to the gradient in velocity
1577     ! of phase M  (shear viscosity in the Mth solids phase)
1578                 MU_s_v(IJK) = MU_sM_sum + MU_sL_ip(IJK,M,M)
1579                 XI_s_v = XI_sM_sum + XI_sL_ip(IJK,M,M)
1580     
1581     ! Bulk viscosity in the Mth solids phase
1582                 LAMBDA_s_v(IJK) = -(2.d0/3.d0)*Mu_s_v(IJK) + XI_s_v
1583     
1584     
1585     ! Find the granular conductivity
1586                 K_s_sum = ZERO
1587     
1588                 K_s_dil = (75.d0/384.d0)*D_PM* RO_S(IJK,M)*&
1589                      DSQRT(PI*Theta_m(IJK,M)/M_PM)
1590     
1591                 IF(SWITCH == ZERO .OR. RO_G(IJK) == ZERO) THEN
1592                    Kth_star = K_s_dil
1593                 ELSEIF(Theta_m(IJK,M)/M_PM < SMALL_NUMBER)THEN
1594                    Kth_star = ZERO
1595     
1596                 ELSEIF(EP_S(IJK,M) <= DIL_EP_s) THEN
1597                    Kth_star = K_s_dil*EP_s(IJK,M)*G_0(IJK,M,M)/ &
1598                        (SUM_EpsGo+ 1.2d0*DgA_sM*K_s_dil &
1599                        / (RO_S(IJK,M)**2 *(Theta_m(IJK,M)/M_PM)))
1600                 ELSE
1601                    Kth_star = K_s_dil*EP_S(IJK,M)*G_0(IJK,M,M)/ &
1602                        (SUM_EpsGo+ 1.2d0*F_gs(IJK,M)*K_s_dil &
1603                        / (RO_S(IJK,M)**2 *EP_s(IJK,M)*(Theta_m(IJK,M)/M_PM)))
1604                 ENDIF
1605     
1606     ! Kth doesn't include the mass.
1607                 K_s_MM = (Kth_star/(M_PM*G_0(IJK,M,M)))*&
1608                      (1.d0+(3.d0/5.d0)*(1.d0+C_E)*(1.d0+C_E)*SUM_EpsGo)**2
1609     
1610                 DO L = 1, SMAX
1611                    D_PL = D_P(IJK,L)
1612                    M_PL = (PI/6.d0)*D_PL**3 *RO_S(IJK,L)
1613                    MPSUM = M_PM + M_PL
1614                    DPSUMo2 = (D_PM+D_PL)/2.d0
1615                    NU_PL = ROP_S(IJK,L)/M_PL
1616     
1617                    IF ( L .eq. M) THEN
1618     
1619     ! solids phase 'conductivity' associated with the
1620     ! difference in velocity. again these terms cancel when
1621     ! added together so do not explicity include them
1622     ! in calculations
1623                       Kvel_s_ip(IJK,M,L) = ZERO
1624                       !     K_common_term*NU_PM*NU_PL*&
1625                       !     (3.d0*PI/10.d0)*R0p_lm*Theta_m(IJK,M)
1626     
1627     ! solids phase 'conductivity' associated with the
1628     ! difference in the gradient in number densities.
1629     ! again these terms cancel so do not explicity include
1630     ! them in calculations
1631                       Knu_sL_ip(IJK,M,L) = ZERO
1632                       !     K_common_term*NU_PM*&
1633                       !     (PI*DPSUMo2/6.d0)*R1p_lm*(Theta_m(IJK,M)**(3./2.))
1634     
1635                       Knu_sM_ip(IJK,M,L) = ZERO
1636                       !     K_common_term*NU_PL*&
1637                       !     (PI*DPSUMo2/6.d0)*R1p_lm*(Theta_m(IJK,M)**(3./2.))
1638     
1639                       K_s_sum = K_s_sum + K_s_MM
1640     
1641                    ELSE
1642                       Ap_lm = (M_PM*Theta_m(IJK,L)+M_PL*Theta_m(IJK,M))/2.d0
1643                       Bp_lm = (M_PM*M_PL*(Theta_m(IJK,L)-&
1644                           Theta_m(IJK,M) ))/(2.d0*MPSUM)
1645                       Dp_lm = (M_PL*M_PM*(M_PM*Theta_m(IJK,M)+&
1646                           M_PL*Theta_m(IJK,L) ))/(2.d0*MPSUM*MPSUM)
1647                       R0p_lm = (1.d0/(Ap_lm**1.5 * Dp_lm**2.5))+&
1648                           ((15.d0*Bp_lm*Bp_lm)/(2.d0* Ap_lm**2.5 * Dp_lm**3.5))+&
1649                           ((175.d0*(Bp_lm**4))/(8.d0*Ap_lm**3.5 * Dp_lm**4.5))
1650     
1651                       R1p_lm = (1.d0/((Ap_lm**1.5)*(Dp_lm**3)))+ &
1652                           ((9.d0*Bp_lm*Bp_lm)/(Ap_lm**2.5 * Dp_lm**4))+&
1653                           ((30.d0*Bp_lm**4)/(2.d0*Ap_lm**3.5 * Dp_lm**5))
1654     
1655                       R5p_lm = (1.d0/(Ap_lm**2.5 * Dp_lm**3 ) )+ &
1656                           ((5.d0*Bp_lm*Bp_lm)/(Ap_lm**3.5 * Dp_lm**4))+&
1657                           ((14.d0*Bp_lm**4)/(Ap_lm**4.5 * Dp_lm**5))
1658     
1659                       R6p_lm = (1.d0/(Ap_lm**3.5 * Dp_lm**3))+ &
1660                           ((7.d0*Bp_lm*Bp_lm)/(Ap_lm**4.5 * Dp_lm**4))+&
1661                           ((126.d0*Bp_lm**4)/(5.d0*Ap_lm**5.5 * Dp_lm**5))
1662     
1663                       R7p_lm = (3.d0/(2.d0*Ap_lm**2.5 * Dp_lm**4))+ &
1664                           ((10.d0*Bp_lm*Bp_lm)/(Ap_lm**3.5 * Dp_lm**5))+&
1665                           ((35.d0*Bp_lm**4)/(Ap_lm**4.5 * Dp_lm**6))
1666     
1667                       R8p_lm = (1.d0/(2.d0*Ap_lm**1.5 * Dp_lm**4))+ &
1668                           ((6.d0*Bp_lm*Bp_lm)/(Ap_lm**2.5 * Dp_lm**5))+&
1669                           ((25.d0*Bp_lm**4)/(Ap_lm**3.5 * Dp_lm**6))
1670     
1671                       R9p_lm = (1.d0/(Ap_lm**2.5 * Dp_lm**3))+ &
1672                           ((15.d0*Bp_lm*Bp_lm)/(Ap_lm**3.5 * Dp_lm**4))+&
1673                           ((70.d0*Bp_lm**4)/(Ap_lm**4.5 * Dp_lm**5))
1674     
1675                       K_common_term = DPSUMo2**3 * M_PL*M_PM/(2.d0*MPSUM)*&
1676                            (1.d0+C_E)*G_0(IJK,M,L) * (M_PM*M_PL)**1.5
1677     
1678     ! solids phase 'conductivity' associated with the
1679     ! gradient in granular temperature of species M
1680                       K_s_LM = - K_common_term*NU_PM*NU_PL*(&
1681                        ((DPSUMo2*DSQRT(PI)/16.d0)*(3.d0/2.d0)*Bp_lm*R5p_lm)+&
1682                        ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*PI/6.d0)*&
1683                        (3.d0/2.d0)*R1p_lm)-(&
1684                        ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PM/8.d0)*Bp_lm*R6p_lm)+&
1685                        ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1686                        8.d0)*M_PM*R9p_lm)+&
1687                        ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PL*M_PM/(MPSUM*MPSUM))*&
1688                        M_PL*Bp_lm*R7p_lm)+&
1689                        ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1690                        2.d0)*(M_PM/MPSUM)**2 * M_PL*R8p_lm)+&
1691                        ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PM*M_PL/(2.d0*MPSUM))*&
1692                        R9p_lm)-&
1693                        ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*DPSUMo2*DSQRT(PI)*&
1694                        (M_PM*M_PL/MPSUM)*Bp_lm*R7p_lm) )*Theta_M(IJK,L) )*&
1695                        (Theta_M(IJK,M)**2 * Theta_M(IJK,L)**3)
1696     
1697     ! solids phase 'conductivity' associated with the
1698     ! gradient in granular temperature of species L
1699                         !Kth_sL_ip(IJK,M,L) = K_common_term*NU_PM*NU_PL*(&
1700                       Kth_sl_iptmp = K_common_term*NU_PM*NU_PL*(&
1701                         (-(DPSUMo2*DSQRT(PI)/16.d0)*(3.d0/2.d0)*Bp_lm*R5p_lm)-&
1702                         ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*PI/6.d0)*&
1703                         (3.d0/2.d0)*R1p_lm)+(&
1704                         ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PL/8.d0)*Bp_lm*R6p_lm)+&
1705                         ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1706                         8.d0)*M_PL*R9p_lm)+&
1707                         ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PL*M_PM/(MPSUM*MPSUM))*&
1708                         M_PM*Bp_lm*R7p_lm)+&
1709                         ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1710                         2.d0)*(M_PM/MPSUM)**2 * M_PM*R8p_lm)+&
1711                         ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PM*M_PL/(2.d0*MPSUM))*&
1712                         R9p_lm)+&
1713                         ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*DPSUMo2*DSQRT(PI)*&
1714                         (M_PM*M_PL/MPSUM)*Bp_lm*R7p_lm) )*Theta_M(IJK,M) )*&
1715                         (Theta_M(IJK,L)**2 * Theta_M(IJK,M)**3)
1716     
1717                       Kth_sL_ip(IJK,M,L) = (ONE-UR_Kth_sml)*Kth_sL_ip(IJK,M,L) +&
1718                         UR_Kth_sml * Kth_sl_iptmp
1719     
1720     ! solids phase 'conductivity' associated with the
1721     ! difference in velocity
1722                       Kvel_s_ip(IJK,M,L) = K_common_term*NU_PM*NU_PL*&
1723                         (M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(3.d0*PI/10.d0)*&
1724                         R0p_lm * (Theta_m(IJK,M)*Theta_m(IJK,L))**2.5
1725     
1726     ! solids phase 'conductivity' associated with the
1727     ! difference in the gradient in number densities
1728                       Knu_sL_ip(IJK,M,L) = K_common_term*(&
1729                         ((DPSUMo2*DSQRT(PI)/16.d0)*Bp_lm*R5p_lm)+&
1730                         ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*PI/6.d0)*&
1731                         R1p_lm) )* (Theta_m(IJK,M)*Theta_m(IJK,L))**3
1732     
1733     ! to avoid recomputing Knu_sL_ip, sof.
1734                       Knu_sM_ip(IJK,M,L) = NU_PL * Knu_sL_ip(IJK,M,L)
1735                       Knu_sL_ip(IJK,M,L) = NU_PM * Knu_sL_ip(IJK,M,L)
1736                       K_s_sum = K_s_sum + K_s_LM
1737                    ENDIF
1738                 ENDDO
1739     
1740     ! granular conductivity in Mth solids phase
1741                 Kth_s(IJK,M) = K_s_sum
1742     
1743              ENDIF   ! Fluid_at
1744           ENDDO   ! IJK loop
1745           RETURN
1746           END SUBROUTINE GT_PDE_IA
1747     
1748     
1749     
1750     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1751     !                                                                      C
1752     !  Subroutine: FRICTION_SCHAEFFER                                      C
1753     !  Purpose: Calculate frictional-flow stress terms                     C
1754     !                                                                      C
1755     !  Author: M. Syamlal                               Date: 10-FEB-93    C
1756     !                                                                      C
1757     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1758           Subroutine friction_schaeffer(M)
1759     
1760     ! Modules
1761     !---------------------------------------------------------------------//
1762     
1763           USE param1, only: zero, small_number
1764     
1765           USE fldvar, only: ep_g, ep_s
1766           USE fldvar, only: p_star, theta_m
1767     
1768           USE physprop, only: smax, close_packed
1769     
1770           USE constant, only: to_SI
1771           USE constant, only: sin2_phi
1772     
1773           use run, only: granular_energy
1774     
1775           USE visc_s, only: ep_g_blend_end
1776           USE visc_s, only: mu_s_p, lambda_s_p
1777           USE visc_s, only: i2_devD_s
1778     
1779           USE compar, only: ijkstart3, ijkend3
1780           USE functions, only: fluid_at
1781           IMPLICIT NONE
1782     
1783     ! Dummy arguments
1784     !---------------------------------------------------------------------//
1785     ! solids phase index
1786           INTEGER, INTENT(IN) :: M
1787     
1788     ! Local parameters
1789     !---------------------------------------------------------------------//
1790     ! maximum value of solids viscosity in poise
1791           DOUBLE PRECISION :: MAX_MU_s
1792           PARAMETER (MAX_MU_s = 1000.D0)
1793     
1794     ! Local variables
1795     !---------------------------------------------------------------------//
1796     ! cell index
1797           INTEGER :: IJK
1798     ! solids phase index
1799           INTEGER :: MM
1800     ! sum of all solids volume fractions
1801           DOUBLE PRECISION :: SUM_EPS_CP
1802     ! factor in frictional-flow stress terms
1803           DOUBLE PRECISION :: qxP_s
1804     !---------------------------------------------------------------------//
1805     
1806           DO IJK = ijkstart3, ijkend3
1807              IF ( FLUID_AT(IJK) ) THEN
1808     
1809                 IF(EP_g(IJK) .LT. EP_g_blend_end(IJK)) THEN
1810     ! Tardos, PT, (1997), pp. 61-74 explains in his equation (3) that
1811     ! solids normal and shear frictional stresses have to be treated
1812     ! consistently. Therefore, add closed_packed check for consistency
1813     ! with the treatment of the normal frictional force (see for example
1814     ! source_v_s.f). sof May 24 2005.
1815                    SUM_EPS_CP=0.0
1816                    DO MM=1,SMAX
1817                       IF (CLOSE_PACKED(MM)) SUM_EPS_CP = SUM_EPS_CP + &
1818                          EP_S(IJK,MM)
1819                    ENDDO
1820     
1821     ! plastic pressure (p_star) is calculated seperately to ensure
1822     ! its value is up-to-date with latest value of ep_g
1823     !               P_star(IJK) = Neg_H(EP_g(IJK),EP_star_array(IJK))
1824     
1825     !-----------------------------------------------------------------------
1826     ! Gray and Stiles (1988) frictional-flow stress tensor
1827     !              IF(Sin2_Phi .GT. SMALL_NUMBER) THEN
1828     !                 qxP_s = SQRT( (4. * Sin2_Phi) * I2_devD_s(IJK) + &
1829     !                    trD_s_C(IJK,M) * trD_s_C(IJK,M))
1830     !                 MU_s_p(IJK) = P_star(IJK) * Sin2_Phi/ &
1831     !                    (qxP_s + SMALL_NUMBER)*(EP_S(IJK,M)/SUM_EPS_CP)
1832     !                 MU_s_p(IJK) = MIN(MU_s_p(IJK), to_SI*MAX_MU_s)
1833     !                 LAMBDA_s_p(IJK) = P_star(IJK) * F_Phi /&
1834     !                    (qxP_s + SMALL_NUMBER)*(EP_S(IJK,M)/SUM_EPS_CP)
1835     !                 LAMBDA_s_p(IJK) = MIN(LAMBDA_s(IJK, M), to_SI*MAX_MU_s)
1836     !              ENDIF
1837     !-----------------------------------------------------------------------
1838     
1839     ! Schaeffer (1987)
1840                    qxP_s = SQRT( (4.D0 * Sin2_Phi) * I2_devD_s(IJK))
1841     ! multiply by factor ep_s/sum_eps_cp for consistency with solids
1842     ! pressure treatment
1843                    MU_s_p(IJK) = P_star(IJK) * Sin2_Phi/&
1844                       (qxP_s + SMALL_NUMBER) * (EP_S(IJK,M)/SUM_EPS_CP)
1845                    MU_s_p(IJK) = MIN(MU_s_p(IJK), to_SI*MAX_MU_s)
1846     
1847                    LAMBDA_s_p(IJK) = ZERO
1848     
1849     ! when solving for the granular energy equation (PDE) setting
1850     ! theta = 0 is done in solve_granular_energy.f to avoid
1851     ! convergence problems. (sof)
1852                    IF(.NOT.GRANULAR_ENERGY) THETA_m(IJK, M) = ZERO
1853                 ENDIF
1854              ENDIF   ! end if (fluid_at(ijk)
1855     
1856           ENDDO   ! IJK loop
1857           RETURN
1858           END SUBROUTINE FRICTION_SCHAEFFER
1859     
1860     
1861     
1862     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1863     !                                                                      C
1864     !  Subroutine: FRICTION_PRINCETON                                      C
1865     !  Purpose: Calculate frictional contributions to granular stress      C
1866     !  terms based on model by Srivastava & Sundaresan (2003)              C
1867     !                                                                      C
1868     !  Author: Anuj Srivastava, Princeton University    Date: 20-APR-98    C
1869     !                                                                      C
1870     !  Literature/Document References:                                     C
1871     !  Srivastava, A., and Sundaresan, S., Analysis of a frictional-       C
1872     !     kinetic model for gas-particle flow, Powder Technology,          C
1873     !     2003, 129, 72-85.                                                C
1874     !  Benyahia, S., Validation study of two continuum granular            C
1875     !     frictional flow theories, Industrial & Engineering Chemistry     C
1876     !     Research, 2008, 47, 8926-8932.                                   C
1877     !                                                                      C
1878     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1879           Subroutine Friction_princeton(M)
1880     
1881     ! Modules
1882     !---------------------------------------------------------------------//
1883           USE constant, only: alpha
1884           USE constant, only: eta
1885           USE constant, only: sqrt_pi
1886           USE constant, only: to_si
1887           USE constant, only: sin_phi
1888           USE constant, only: eps_f_min
1889           USE constant, only: fr, n_pc, d_pc, n_pf, delta
1890     
1891           USE fldvar, only: ep_g, ep_s
1892           USE fldvar, only: d_p, ro_s
1893           USE fldvar, only: theta_m
1894           USE fldvar, only: p_s_f
1895     
1896           USE param1, only: zero, one, small_number
1897     
1898           USE physprop, only: smax, close_packed
1899     
1900           USE rdf, only: g_0
1901     
1902           USE run, only: kt_type_enum, ghd_2007
1903           USE run, only: savage
1904           USE trace, only: trd_s2, trd_s_c
1905           USE visc_s, only: mu_s_f, lambda_s_f
1906           USE visc_s, only: mu_s_v, mu_b_v
1907           USE visc_s, only: ep_star_array
1908     
1909           USE compar, only: ijkstart3, ijkend3
1910           USE functions, only: fluid_at
1911           IMPLICIT NONE
1912     
1913     ! Dummy arguments
1914     !---------------------------------------------------------------------//
1915     ! solids phase index
1916           INTEGER, INTENT(IN) :: M
1917     
1918     ! Local variables
1919     !---------------------------------------------------------------------//
1920     ! cell index
1921           INTEGER :: IJK
1922     ! solids phase index
1923           INTEGER :: MM
1924     ! used to compute frictional terms
1925           DOUBLE PRECISION :: Chi, Pc, Mu_zeta,PfoPc, N_Pff
1926           DOUBLE PRECISION :: ZETA
1927     ! sum of all solids volume fractions
1928           DOUBLE PRECISION :: SUM_EPS_CP
1929     !     parameters in pressure linearization; simple averaged Dp
1930           DOUBLE PRECISION :: dpc_dphi, dp_avg
1931     !---------------------------------------------------------------------//
1932     
1933           DO IJK = ijkstart3, ijkend3
1934              IF ( FLUID_AT(IJK) ) THEN
1935     
1936                 IF (EP_g(IJK) .LT. (ONE-eps_f_min)) THEN
1937     
1938     ! This close_packed check was added for consistency with the Schaeffer
1939     ! model implementation.
1940                    SUM_EPS_CP = ZERO
1941                    dp_avg = ZERO
1942                    DO MM=1,SMAX
1943                       dp_avg = dp_avg + D_p(IJK,MM)
1944                       IF (CLOSE_PACKED(MM)) SUM_EPS_CP = SUM_EPS_CP + &
1945                          EP_S(IJK,MM)
1946                    ENDDO
1947                    dp_avg = dp_avg/DBLE(SMAX)
1948     
1949                    IF (SAVAGE.EQ.1) THEN !form of Savage (not to be used with GHD theory)
1950                       Mu_zeta = &
1951                          ((2d0+ALPHA)/3d0)*((Mu_s_v(IJK)/(Eta*(2d0-Eta)*&
1952                          G_0(IJK,M,M)))*(1d0+1.6d0*Eta*EP_s(IJK,M)*&
1953                          G_0(IJK,M,M))*(1d0+1.6d0*Eta*(3d0*Eta-2d0)*&
1954                          EP_s(IJK,M)*G_0(IJK,M,M))+(0.6d0*Mu_b_v(IJK)*Eta))
1955                       ZETA = &
1956                          ((48d0*Eta*(1d0-Eta)*RO_S(IJK,M)*EP_s(IJK,M)*&
1957                          EP_s(IJK,M)*G_0(IJK,M,M)*&
1958                          (Theta_m(IJK,M)**1.5d0))/&
1959                          (SQRT_Pi*D_p(IJK,M)*2d0*Mu_zeta))**0.5d0
1960     
1961                    ELSEIF (SAVAGE.EQ.0) THEN !S:S form
1962                       ZETA = (SMALL_NUMBER +&
1963                          trD_s2(IJK,M) - ((trD_s_C(IJK,M)*&
1964                          trD_s_C(IJK,M))/3.d0))**0.5d0
1965     
1966                    ELSE          !combined form
1967                       IF(KT_TYPE_ENUM == GHD_2007) THEN
1968                         ZETA = ((Theta_m(IJK,M)/dp_avg**2) +&
1969                            (trD_s2(IJK,M) - ((trD_s_C(IJK,M)*&
1970                            trD_s_C(IJK,M))/3.d0)))**0.5d0
1971                       ELSE
1972                         ZETA = ((Theta_m(IJK,M)/D_p(IJK,M)**2) +&
1973                            (trD_s2(IJK,M) - ((trD_s_C(IJK,M)*&
1974                            trD_s_C(IJK,M))/3.d0)))**0.5d0
1975                       ENDIF
1976                    ENDIF
1977     
1978     
1979                    IF ((ONE-EP_G(IJK)) .GT. ((ONE-ep_star_array(ijk))-delta)) THEN
1980     ! Linearized form of Pc; this is more stable and provides continuous function.
1981                       dpc_dphi = (to_SI*Fr)*((delta**5)*(2d0*(ONE-&
1982                           ep_star_array(IJK)-delta) - 2d0*eps_f_min)+&
1983                           ((ONE-ep_star_array(ijk)-delta)-eps_f_min)*&
1984                           (5*delta**4))/(delta**10)
1985     
1986                       Pc = (to_SI*Fr)*(((ONE-ep_star_array(IJK)-delta) -&
1987                          EPS_f_min)**N_Pc)/(delta**D_Pc)
1988                       Pc = Pc + dpc_dphi*((ONE-EP_G(IJK))+delta-(ONE-&
1989                          ep_star_array(IJK)))
1990     
1991                     ! Pc = 1d25*(((ONE-EP_G(IJK))- (ONE-ep_star_array(ijk)))**10d0) ! old commented Pc
1992                    ELSE
1993                       Pc = (to_SI*Fr)*(((ONE-EP_G(IJK)) - EPS_f_min)**N_Pc) / &
1994                          (((ONE-ep_star_array(ijk)) - (ONE-EP_G(IJK)) +&
1995                          SMALL_NUMBER)**D_Pc)
1996                    ENDIF
1997     
1998                    IF (trD_s_C(IJK,M) .GE. ZERO) THEN
1999                       N_Pff = DSQRT(3d0)/(2d0*Sin_Phi) !dilatation
2000                    ELSE
2001                       N_Pff = N_Pf !compaction
2002                    ENDIF
2003     
2004                    IF ((trD_s_C(IJK,M)/(ZETA*N_Pff*DSQRT(2d0)&
2005                         *Sin_Phi)) .GT. 1d0) THEN
2006                       P_s_f(IJK) =ZERO
2007                       PfoPc = ZERO
2008                    ELSEIF(trD_s_C(IJK,M) == ZERO) THEN
2009                       P_s_f(IJK) = Pc
2010                       PfoPc = ONE
2011                    ELSE
2012                       P_s_f(IJK) = Pc*(1d0 - (trD_s_C(IJK,M)/(ZETA&
2013                          *N_Pff*DSQRT(2d0)*Sin_Phi)))**(N_Pff-1d0)
2014                       PfoPc = (1d0 - (trD_s_C(IJK,M)/(ZETA&
2015                          *N_Pff*DSQRT(2d0)*Sin_Phi)))**(N_Pff-1d0)
2016                    ENDIF
2017     
2018                    Chi = DSQRT(2d0)*P_s_f(IJK)*Sin_Phi*(N_Pff - (N_Pff-1d0)*&
2019                       (PfoPc)**(1d0/(N_Pff-1d0)))
2020     
2021                    IF (Chi < ZERO) THEN
2022                       P_s_f(IJK) = Pc*((N_Pff/(N_Pff-1d0))**(N_Pff-1d0))
2023                       Chi = ZERO
2024                    ENDIF
2025     
2026                    Mu_s_f(IJK) = Chi/(2d0*ZETA)
2027                    Lambda_s_f(IJK) = -2d0*Mu_s_f(IJK)/3d0
2028     
2029     ! Modification of the stresses when more than one solids phase is used:
2030     ! This is NOT done when mixture mom. eq. is solved (i.e. for GHD theory)
2031                    IF(KT_TYPE_ENUM /= GHD_2007) THEN
2032                       P_s_f(IJK) = P_s_f(IJK) * (EP_S(IJK,M)/SUM_EPS_CP)
2033                       Mu_s_f(IJK) = Mu_s_f(IJK) * (EP_S(IJK,M)/SUM_EPS_CP)
2034                       Lambda_s_f(IJK) = Lambda_s_f(IJK) * (EP_S(IJK,M)/SUM_EPS_CP)
2035                    ENDIF
2036     
2037                 ENDIF   ! end if ep_g < 1-eps_f_min
2038              ENDIF   ! Fluid_at
2039           ENDDO   ! IJK loop
2040           RETURN
2041           END SUBROUTINE FRICTION_PRINCETON
2042     
2043     
2044     
2045     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2046     !                                                                      C
2047     !  Subroutine: SUBGRID_STRESS_IGCI                                     C
2048     !  Purpose: Calculate solids viscosity and pressure using subgrid      C
2049     !  model                                                               C
2050     !                                                                      C
2051     !  Author: Sebastien Dartevelle, LANL, May 2013                        C
2052     !                                                                      C
2053     !  Literature/Document References:                                     C
2054     !  Igci, Y., Pannala, S., Benyahia, S., & Sundaresan S.,               C
2055     !     Validation studies on filtered model equations for gas-          C
2056     !     particle flows in risers, Industrial & Engineering Chemistry     C
2057     !     Research, 2012, 51(4), 2094-2103                                 C
2058     !                                                                      C
2059     !  Comments:                                                           C
2060     !  Still needs to be reviewed for accuracy with source material        C
2061     !                                                                      C
2062     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2063           Subroutine subgrid_stress_igci(M)
2064     
2065     ! Modules
2066     !---------------------------------------------------------------------//
2067           USE param1, only: zero, one, small_number
2068     
2069           USE fldvar, only: ep_s
2070           USE fldvar, only: d_p, ro_s
2071           USE fldvar, only: ro_g
2072           USE fldvar, only: p_s_v
2073     ! this shouldn't be necessary
2074           use fldvar, only: theta_m
2075     
2076           USE visc_s, only: mu_s_v, lambda_s_v
2077     
2078           USE physprop, only: mu_g
2079     
2080           USE run, only: filter_size_ratio
2081           USE run, only: subgrid_wall
2082     
2083           USE constant, only: gravity
2084     
2085           use geometry, only: vol, axy, do_k
2086           USE compar, only: ijkstart3, ijkend3
2087           USE functions, only: fluid_at
2088           IMPLICIT NONE
2089     
2090     ! Dummy arguments
2091     !---------------------------------------------------------------------//
2092     ! solids phase index
2093           INTEGER, INTENT(IN) :: M
2094     
2095     ! Local variables
2096     !---------------------------------------------------------------------//
2097     ! cell indices
2098           INTEGER :: IJK
2099     ! Igci models
2100           DOUBLE PRECISION :: pressurefac, ps_kinetic, ps_total, ps_sub
2101           DOUBLE PRECISION :: viscosityfac, extra_visfactor, mu_sub
2102           DOUBLE PRECISION :: mu_kinetic, mu_total
2103     ! factors to correct subgrid effects from the wall
2104           DOUBLE PRECISION :: wfactor_Ps, wfactor_mus
2105     ! the inverse Froude number, or dimensionless Filtersize
2106           DOUBLE PRECISION :: Inv_Froude
2107     ! one particle terminal settling velocity
2108           DOUBLE PRECISION :: vt
2109     ! the filter size which is a function of each grid cell volume
2110           DOUBLE PRECISION :: filtersize
2111     !---------------------------------------------------------------------//
2112     
2113           DO IJK = ijkstart3, ijkend3
2114              IF ( FLUID_AT(IJK) ) THEN
2115     
2116     ! initialize
2117                 wfactor_Ps = ONE   ! for P_S
2118                 wfactor_mus = ONE   ! for MU_s
2119     
2120     ! particle terminal settling velocity: vt = g*d^2*(Rho_s - Rho_g) / 18 * Mu_g
2121                 vt = GRAVITY*D_p(IJK,M)*D_p(IJK,M)*&
2122                      (RO_S(IJK,M) - RO_g(IJK)) / (18.0d0*MU_G(IJK))
2123     
2124     ! FilterSIZE calculation for each specific gridcell volume
2125                 IF(DO_K) THEN
2126                    filtersize = filter_size_ratio * (VOL(IJK)**(ONE/3.0d0))
2127                 ELSE
2128                    filtersize = filter_size_ratio * DSQRT(AXY(IJK))
2129                 ENDIF
2130     
2131     
2132     ! dimensionless inverse of Froude number
2133                 Inv_Froude =  filtersize * GRAVITY / vt**2
2134     
2135     ! various factor needed:
2136                 pressurefac = 0.48d0*(Inv_Froude**0.86)*&
2137                    (ONE-EXP(-Inv_Froude/1.4))
2138                 viscosityfac = 0.37d0*(Inv_Froude**1.22)
2139                 Extra_visfactor = ONE/(0.28d0*(Inv_Froude**0.43)+ONE)
2140     
2141                 IF (EP_s(IJK,M) .LE. 0.0131) THEN
2142                    Ps_kinetic = -10.4d0*(EP_s(IJK,m)**2)+0.31d0*EP_s(IJK,m)
2143                 ELSEIF (EP_s(IJK,M) .LE. 0.290) THEN
2144                    Ps_kinetic = -0.185d0*(EP_s(IJK,m)**3)+&
2145                       0.066d0*(EP_s(IJK,m)**2)-0.000183d0*EP_s(IJK,m)+&
2146                       0.00232d0
2147                 ELSEIF (EP_s(IJK,M) .LE. 0.595) THEN
2148                    Ps_kinetic = -0.00978d0*EP_s(IJK,m)+0.00615d0
2149                 ELSE
2150                    Ps_kinetic = -6.62d0*(EP_s(IJK,m)**3)+&
2151                       49.5d0*(EP_s(IJK,m)**2)-50.3d0*EP_s(IJK,m)+13.8d0
2152                 ENDIF
2153     
2154                 Ps_sub = pressurefac*(EP_s(IJK,M)-0.59d0)*&
2155                    (-1.69d0*EP_s(IJK,M)-4.61d0*(EP_s(IJK,M)**2)+&
2156                    11.d0*(EP_s(IJK,M)**3))
2157     
2158                 IF (Ps_sub .GE. ZERO) THEN
2159                    Ps_total=Ps_kinetic+Ps_sub
2160                 ELSE
2161                    Ps_total=Ps_kinetic
2162                 ENDIF
2163     
2164                 IF (EP_s(IJK,M) .LE. 0.02) THEN
2165                    Mu_kinetic = 1720.d0*(EP_s(IJK,m)**4)-&
2166                       215.d0*(EP_s(IJK,m)**3) + 9.81d0*(EP_s(IJK,m)**2)-&
2167                       0.207d0*EP_s(IJK,m)+0.00254d0
2168                 ELSEIF (EP_s(IJK,M) .LE. 0.2) THEN
2169                    Mu_kinetic = 2.72d0*(EP_s(IJK,m)**4)-&
2170                       1.55d0*(EP_s(IJK,m)**3)+0.329d0*(EP_s(IJK,m)**2)-&
2171                       0.0296d0*EP_s(IJK,m)+0.00136d0
2172                 ELSEIF (EP_s(IJK,M) .LE. 0.6095) THEN
2173                    Mu_kinetic = -0.0128d0*(EP_s(IJK,m)**3)+&
2174                       0.0107d0*(EP_s(IJK,m)**2)-0.0005d0*EP_s(IJK,m)+&
2175                       0.000335d0
2176                 ELSE
2177                    Mu_kinetic = 23.6d0*(EP_s(IJK,m)**2)-&
2178                       28.0d0*EP_s(IJK,m)+8.30d0
2179                 ENDIF
2180     
2181                 Mu_sub = Extra_visfactor*viscosityfac*&
2182                    (EP_s(IJK,M)-0.59d0)*(-1.22d0*EP_s(IJK,M)-&
2183                    0.7d0*(EP_s(IJK,M)**2)-2.d0*(EP_s(IJK,M)**3))
2184     
2185                 IF (Mu_sub .GE. ZERO) THEN
2186                    Mu_total = Mu_kinetic+Mu_sub
2187                 ELSE
2188                    Mu_total = Mu_kinetic
2189                 ENDIF
2190     
2191                 IF (SUBGRID_WALL) THEN
2192                    CALL SUBGRID_STRESS_WALL(wfactor_Ps,wfactor_Mus,vt,IJK)
2193                 ENDIF
2194     
2195     ! pressure
2196                 P_s_v(IJK) = Ps_total * wfactor_Ps * (vt**2) * &
2197                    RO_S(IJK,M)
2198     
2199     ! shear viscosity
2200                 Mu_s_v(IJK) = Mu_total * wfactor_mus * (vt**3) * &
2201                    RO_S(IJK,M)/GRAVITY
2202     
2203     ! set an arbitrary value in case value gets negative (this should
2204     ! not happen unless filtersize becomes unrelastic w.r.t. gridsize)
2205                 IF (P_s_v(IJK) .LE. SMALL_NUMBER) P_s_v(IJK) = SMALL_NUMBER
2206                 IF (Mu_s_v(IJK) .LE. SMALL_NUMBER) Mu_s_v(IJK)= SMALL_NUMBER
2207     
2208     ! Solid second viscosity, assuming the bulk viscosity is ZERO
2209                 lambda_s_v(IJK) = (-2.0d0/3.0d0)*Mu_s_v(IJK)
2210     
2211     ! granular temperature is zeroed in all LES/Subgrid model
2212     ! why is this necessary? theta_M shouldn't be used anyway!
2213                 THETA_m(IJK, M) = ZERO
2214     
2215              ENDIF   ! endif (fluid_at(IJK))
2216           ENDDO   ! outer IJK loop
2217           RETURN
2218           END SUBROUTINE SUBGRID_STRESS_IGCI
2219     
2220     
2221     
2222     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2223     !                                                                      C
2224     !  Subroutine: SUBGRID_STRESS_MILIOLI                                  C
2225     !  Purpose: Calculate solids viscosity and pressure using subgrid      C
2226     !  model                                                               C
2227     !                                                                      C
2228     !  Author: Sebastien Dartevelle, LANL, May 2013                        C
2229     !                                                                      C
2230     !  Literature/Document References:                                     C
2231     !  Milioli, C. C., et al., Filtered two-fluid models of fluidized      C
2232     !     gas-particle flows: new constitutive relations, AICHE J,         C
2233     !     doi: 10.1002/aic.14130                                           C
2234     !                                                                      C
2235     !  Comments:                                                           C
2236     !  Still needs to be reviewed for accuracy with source material        C
2237     !                                                                      C
2238     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2239           Subroutine subgrid_stress_MILIOLI(M)
2240     
2241     ! Modules
2242     !---------------------------------------------------------------------//
2243           USE param1, only: zero, one, small_number
2244     
2245           USE fldvar, only: ep_g, ep_s
2246           USE fldvar, only: d_p, ro_s
2247           USE fldvar, only: ro_g
2248           USE fldvar, only: p_s_v
2249     ! this shouldn't be necessary
2250           use fldvar, only: theta_m
2251     
2252           USE visc_s, only: mu_s_v, lambda_s_v
2253           USE visc_s, only: i2_devD_s
2254     
2255           USE physprop, only: mu_g
2256     
2257           USE run, only: filter_size_ratio
2258           USE run, only: subgrid_wall
2259     
2260           USE constant, only: gravity
2261     
2262           use geometry, only: vol, axy, do_k
2263           USE compar, only: ijkstart3, ijkend3
2264           USE functions, only: fluid_at
2265           IMPLICIT NONE
2266     
2267     ! Dummy arguments
2268     !---------------------------------------------------------------------//
2269     ! solids phase index
2270           INTEGER, INTENT(IN) :: M
2271     
2272     ! Local variables
2273     !---------------------------------------------------------------------//
2274     ! cell indices
2275           INTEGER :: IJK
2276     ! Milioli model
2277           DOUBLE PRECISION :: cvisc_pot, cvisc_num, cvisc_den, Cvisc
2278           DOUBLE PRECISION :: cpress_pot, cpress_num, cpress_den, Cpress
2279     ! factor functions to correct subgrid effects from the wall
2280           DOUBLE PRECISION :: wfactor_Ps, wfactor_mus
2281     ! the inverse Froude number, or dimensionless Filtersize
2282           DOUBLE PRECISION :: Inv_Froude
2283     ! one particle terminal settling velocity
2284           DOUBLE PRECISION :: vt
2285     ! the filter size which is a function of each grid cell volume
2286           DOUBLE PRECISION :: filtersize
2287     !---------------------------------------------------------------------//
2288     
2289           DO IJK = ijkstart3, ijkend3
2290              IF ( FLUID_AT(IJK) ) THEN
2291     
2292     ! initialize
2293                 wfactor_Ps = ONE   ! for P_S
2294                 wfactor_mus = ONE   ! for MU_s
2295     
2296     ! particle terminal settling velocity: vt = g*d^2*(Rho_s - Rho_g) / 18 * Mu_g
2297                 vt = GRAVITY*D_p(ijk,M)*D_p(ijk,M)*&
2298                      (RO_S(IJK,M)-RO_g(IJK)) / (18.0d0*MU_G(IJK))
2299     
2300     ! FilterSIZE calculation for each specific gridcell volume
2301                 IF(DO_K) THEN
2302                    filtersize = filter_size_ratio * (VOL(IJK)**(ONE/3.0d0))
2303                 ELSE
2304                    filtersize = filter_size_ratio * DSQRT(AXY(IJK))
2305                 ENDIF
2306     
2307     ! dimensionless inverse of Froude number
2308                 Inv_Froude =  filtersize * GRAVITY / vt**2
2309     
2310     ! Cvisc:
2311                 cvisc_pot = (0.59d0-(ONE-EP_g(IJK)))
2312                 cvisc_num = (0.7d0*(ONE-EP_g(IJK))*cvisc_pot)
2313                 cvisc_den = (0.8d0+(17.d0*cvisc_pot*cvisc_pot*cvisc_pot))
2314                 IF ((ONE-EP_g(IJK)) .GE. ZERO .AND. &
2315                     (ONE-EP_g(IJK)) .LE. 0.59) THEN
2316                    cvisc=(cvisc_num/cvisc_den)
2317                 ELSE
2318                    cvisc=ZERO
2319                 ENDIF
2320     
2321     ! aCvisc:
2322     !            IF ((ONE-EP_g(IJK)) .GE. ZERO .AND. &
2323     !                (ONE-EP_g(IJK)) .LE. 0.59) THEN
2324     !               acvisc(IJK) = (0.7d0*(ONE-EP_g(IJK))*(0.59d0-&
2325     !                  (ONE-EP_g(IJK))))/(0.8d0+17.d0*(0.59d0-&
2326     !                  (ONE-EP_g(IJK)))*(0.59d0-(ONE-EP_g(IJK)))*&
2327     !                  (0.59d0-(ONE-EP_g(IJK))))
2328     !            ELSE
2329     !               acvisc(IJK)=ZERO
2330     !            ENDIF
2331     
2332     ! Cpress:
2333                 cpress_pot = (0.59d0-(ONE-EP_g(IJK)))
2334                 cpress_num = (0.4d0*(ONE-EP_g(IJK))*cpress_pot)
2335                 cpress_den = (0.5d0+(13.d0*cpress_pot*cpress_pot*cpress_pot))
2336                 IF ((ONE-EP_g(IJK)) .GE. ZERO .AND. &
2337                     (ONE-EP_g(IJK)) .LE. 0.59) THEN
2338                    cpress = (cpress_num/cpress_den)
2339                 ELSE
2340                    cpress = ZERO
2341                 ENDIF
2342     
2343     ! aCpress
2344     !            IF ((ONE-EP_g(IJK)) .GE. ZERO .AND. &
2345     !                (ONE-EP_g(IJK)) .LE. 0.59) THEN
2346     !                acpress(IJK) = (0.4d0*(ONE-EP_g(IJK))*(0.59d0-&
2347     !                   (ONE-EP_g(IJK))))/(0.5d0+13.d0*(0.59d0-&
2348     !                   (ONE-EP_g(IJK)))*(0.59d0-(ONE-EP_g(IJK)))*&
2349     !                   (0.59d0-(ONE-EP_g(IJK))))
2350     !            ELSE
2351     !            acpress(IJK)=ZERO
2352     !            ENDIF
2353     
2354                 IF (SUBGRID_WALL) THEN
2355                    CALL SUBGRID_STRESS_WALL(wfactor_Ps,wfactor_Mus,vt,IJK)
2356                 ENDIF
2357     
2358     ! solid filtered pressure
2359                 P_s_v(IJK) = RO_S(IJK,M) * Inv_froude**(2/7) * &
2360                    filtersize**2 * DSQRT( I2_devD_s(IJK) )**2 * &
2361                    cpress * wfactor_Ps    !16/7-2=2/7 in [Pa or kg/m.s2]
2362     
2363     ! solids filtered shear viscosity
2364                 Mu_s_v(IJK) = RO_S(IJK,M) * filtersize**2 * &
2365                    DSQRT( I2_devD_s(IJK) ) * cvisc * wfactor_mus  ! [kg/m.s]
2366     
2367     ! set an arbitrary value in case value gets negative (this should
2368     ! not happen unless filtersize becomes unrealistic w.r.t. gridsize)
2369                 IF (P_s_v(IJK) .LE. SMALL_NUMBER) P_s_v(IJK) = SMALL_NUMBER
2370                 IF (Mu_s_v(IJK) .LE. SMALL_NUMBER) Mu_s_v(IJK)= SMALL_NUMBER
2371     
2372     ! Solids second viscosity, assuming the bulk viscosity is ZERO
2373                 lambda_s_v(IJK) = (-2.0d0/3.0d0)*Mu_s_v(IJK)
2374     
2375     ! granular temperature is zeroed in all LES/Subgrid model
2376     ! why is this necessary? theta_M shouldn't be used anyway!
2377                 THETA_m(IJK, M) = ZERO
2378     
2379              ENDIF   ! endif (fluid_at(IJK))
2380           ENDDO   ! outer IJK loop
2381           RETURN
2382           END SUBROUTINE SUBGRID_STRESS_MILIOLI
2383     
2384     
2385     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2386     !                                                                      C
2387     !  Subroutine: SUBGRID_STRESS_WALL                                     C
2388     !  Purpose: Calculate subgrid corrections arising from wall to solids  C
2389     !  viscosity and pressure.                                             C
2390     !                                                                      C
2391     !  Author: Sebastien Dartevelle, LANL, May 2013                        C
2392     !                                                                      C
2393     !  Literature/Document References:                                     C
2394     !  Igci, Y., and Sundaresan, S., Verification of filtered two-         C
2395     !     fluid models for gas-particle flows in risers, AICHE J.,         C
2396     !     2011, 57 (10), 2691-2707.                                        C
2397     !                                                                      C
2398     !  Comments: Currently only valid for free-slip wall but no checks     C
2399     !  are made to ensure user has selected free-slip wall when this       C
2400     !  option is invoked                                                   C
2401     !                                                                      C
2402     !                                                                      C
2403     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2404           Subroutine subgrid_stress_wall(lfactor_ps, lfactor_mus, vt, &
2405                      IJK)
2406     
2407     ! Modules
2408     !---------------------------------------------------------------------//
2409           USE param1, only: ONE
2410           USE constant, only : GRAVITY
2411           USE cutcell, only : DWALL
2412           IMPLICIT NONE
2413     
2414     ! Dummy arguments
2415     !---------------------------------------------------------------------//
2416     ! factor to correct the solids pressure
2417           DOUBLE PRECISION, INTENT(OUT) :: lfactor_ps
2418     ! factor to correct the solids viscosity
2419           DOUBLE PRECISION, INTENT(OUT) :: lfactor_mus
2420     ! one particle terminal settling velocity
2421           DOUBLE PRECISION, INTENT(IN) :: vt
2422     ! current ijk index
2423           INTEGER, INTENT(IN) :: IJK
2424     
2425     ! Local parameters
2426     !---------------------------------------------------------------------//
2427     ! values are only correct for free-slip walls
2428           DOUBLE PRECISION, PARAMETER :: aps=9.14d0, bps=0.345d0,&
2429                                          amus=5.69d0, bmus=0.228d0
2430     
2431     ! Local variables
2432     !---------------------------------------------------------------------//
2433     ! dimensionless distance to the wall
2434           DOUBLE PRECISION :: x_d
2435     !---------------------------------------------------------------------//
2436     
2437     ! initialize
2438           lfactor_ps = ONE
2439           lfactor_mus = ONE
2440     
2441     ! dimensionless distance to the Wall
2442           x_d = DWALL(IJK) * GRAVITY / vt**2
2443     ! wall function for pressure
2444           lfactor_Ps = ONE / ( ONE + aps * (EXP(-bps*x_d)) )
2445     ! wall function for viscosity
2446           lfactor_mus = ONE / ( ONE + amus * (EXP(-bmus*x_d)) )
2447     
2448           RETURN
2449           END SUBROUTINE SUBGRID_STRESS_WALL
2450     
2451     
2452     
2453     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2454     !                                                                      C
2455     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2456           Subroutine add_shear(M)
2457     
2458     ! Modules
2459     !---------------------------------------------------------------------//
2460     ! y-component of solids velocity
2461           USE fldvar, only: v_s
2462     ! shear velocity
2463           USE vshear, only: vsh
2464           USE compar, only: ijkstart3, ijkend3
2465           USE functions, only: fluid_at
2466           IMPLICIT NONE
2467     
2468     ! Dummy arguments
2469     !---------------------------------------------------------------------//
2470     ! solids phase index
2471           INTEGER, INTENT(IN) :: M
2472     
2473     ! Local variables
2474     !---------------------------------------------------------------------//
2475     ! cell index
2476           INTEGER :: IJK
2477     !---------------------------------------------------------------------//
2478     
2479     !!     $omp parallel do private(IJK)
2480           DO IJK= ijkstart3, ijkend3
2481              IF (FLUID_AT(IJK)) THEN
2482                 V_s(ijk,m)=V_s(IJK,m)+VSH(IJK)
2483              ENDIF
2484           ENDDO
2485     
2486           RETURN
2487           END SUBROUTINE ADD_SHEAR
2488     
2489     
2490     
2491     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2492     !                                                                      C
2493     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2494           Subroutine remove_shear(M)
2495     
2496     ! Modules
2497     !---------------------------------------------------------------------//
2498     ! y-component of solids velocity
2499           USE fldvar, only: v_s
2500     ! shear velocity
2501           USE vshear, only: vsh
2502     ! needed for function.inc
2503           USE compar, only: ijkstart3, ijkend3
2504           USE functions, only: fluid_at
2505           IMPLICIT NONE
2506     
2507     ! Dummy arguments
2508     !---------------------------------------------------------------------//
2509     ! solids phase index
2510           INTEGER, INTENT(IN) :: M
2511     
2512     ! Local variables
2513     !---------------------------------------------------------------------//
2514     ! cell index
2515           INTEGER :: IJK
2516     !---------------------------------------------------------------------//
2517     
2518     !!     $omp parallel do private(IJK)
2519           DO IJK= ijkstart3, ijkend3
2520              IF (FLUID_AT(IJK)) THEN
2521                 V_s(IJK,m)=V_s(IJK,m)-VSH(IJK)
2522              ENDIF
2523           ENDDO
2524     
2525           RETURN
2526           END SUBROUTINE REMOVE_SHEAR
2527     
2528     
2529     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2530     !                                                                      C
2531     ! Subroutine: INIT0_MU_S                                               C
2532     !                                                                      C
2533     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2534           Subroutine init0_mu_s (M)
2535     
2536     ! Modules
2537     !---------------------------------------------------------------------//
2538           USE param1, only: zero
2539           USE constant, only: v_ex
2540           USE fldvar, only: p_s_v, p_s_p, p_s_f
2541           USE visc_s, only: mu_s_v, mu_s_p, mu_s_f, mu_b_v
2542           USE visc_s, only: lambda_s_v, lambda_s_p, lambda_s_f
2543           USE visc_s, only: alpha_s
2544     !---------------------------------------------------------------------//
2545           IMPLICIT NONE
2546           INTEGER, INTENT(IN) :: M
2547     !---------------------------------------------------------------------//
2548     
2549           IF (V_EX /= ZERO)  ALPHA_s(:,M) = ZERO
2550     
2551     ! viscous contributions
2552           Mu_s_v(:)     = ZERO
2553           Mu_b_v(:)     = ZERO
2554           LAMBDA_s_v(:) = ZERO
2555           P_s_v(:) = ZERO
2556     ! plastic contributions (local to this routine)
2557           Mu_s_p(:)     = ZERO
2558           LAMBDA_s_p(:) = ZERO   ! never used
2559           P_s_p(:) = ZERO        ! never used
2560     ! frictional contributions
2561           Mu_s_f(:)     = ZERO
2562           LAMBDA_s_f(:) = ZERO
2563           P_s_f(:) = ZERO
2564     
2565           RETURN
2566           END SUBROUTINE init0_mu_s
2567     
2568     
2569     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2570     !                                                                      C
2571     ! Subroutine: INIT1_MU_S                                               C
2572     ! Purpose: Calculate quantities that are functions of velocity only!   C
2573     !  - strain rate tensor: D_s (local)                                   C
2574     !  - trace of strain rate tensor: trd_s_c (global)                     C
2575     !    this quantity seems equivalent to the variable trd_s that is      C
2576     !    calculated by the subroutine calc_trd_s (trd_s). formulation is   C
2577     !    equivalent despite apparent differences:                          C
2578     !       trd (D) = du/dx + dv/dy + 1/x dw/dz + u/x (here)               C
2579     !               = 1/x d(xu)/dx + dv/dy + 1/x dw_dz (calc_trd_s)        C
2580     !  - trace of square of strain rate tensor: trd_s2 (global)            C
2581     !  - second invariant of the deviatoric stess tensor: i2_devD_s        C
2582     !    (global)                                                          C
2583     !                                                                      C
2584     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2585           Subroutine init1_mu_s (M)
2586     
2587     ! Modules
2588     !---------------------------------------------------------------------//
2589           USE compar, only: ijkstart3, ijkend3
2590           USE constant, only: v_ex
2591           USE cutcell, only: cut_cell_at
2592           USE functions, only: fluid_at
2593           USE geometry, only: cylindrical
2594           USE physprop, only: smax
2595           USE param1, only: zero, one, half
2596           USE run, only: kt_type_enum
2597           USE run, only: ia_2005
2598           USE run, only: shear
2599           USE trace, only: trD_s_c
2600           USE trace, only: trD_s2
2601           USE trace, only: trD_s2_ip
2602           USE visc_s, only: I2_devD_s
2603           IMPLICIT NONE
2604     
2605     ! Dummy arguments
2606     !---------------------------------------------------------------------//
2607     ! solids phase index
2608           INTEGER, INTENT(IN) :: M
2609     
2610     ! Local variables
2611     !---------------------------------------------------------------------//
2612     ! solids phase index
2613           INTEGER :: L
2614     ! Indices
2615           INTEGER :: IJK
2616     ! Local DO-LOOP counters and phase index
2617           INTEGER :: I1, I2
2618     ! Strain rate tensor components for mth and lth solids phases
2619           DOUBLE PRECISION :: D_sM(3,3), D_sl(3,3)
2620     ! velocity gradient for mth and lth solids phases
2621           DOUBLE PRECISION :: DelV_sM(3,3), DelV_sL(3,3)
2622     ! shear related reciprocal time scale
2623           DOUBLE PRECISION :: SRT
2624     !---------------------------------------------------------------------//
2625     
2626           IF (SHEAR) THEN
2627              CALL add_shear(M)
2628              IF (KT_TYPE_ENUM == IA_2005) THEN
2629                 DO L = 1,SMAX
2630                    IF (L .NE. M) THEN
2631                       CALL add_shear(L)
2632                    ENDIF
2633                 ENDDO
2634              ENDIF
2635           ENDIF
2636     
2637     !!$omp  parallel do default(shared) &
2638     !!$omp  private( IJK, L, I1, I2, D_sM, D_sL, DelV_sM, DelV_sL)
2639           DO IJK = ijkstart3, ijkend3
2640              IF ( FLUID_AT(IJK) ) THEN
2641     
2642     ! Velocity derivatives (gradient and rate of strain tensor) for Mth
2643     ! solids phase at i, j, k
2644                 CALL CALC_DERIV_VEL_SOLIDS(IJK, M, DelV_sM, D_sM)
2645     
2646     ! Calculate the trace of D_s
2647     ! friction/algebraic various kt
2648                 trD_s_C(IJK,M) = D_sM(1,1) + D_sM(2,2) + D_sM(3,3)
2649     
2650     ! Calculate trace of the square of D_s
2651     ! friction/algebraic various kt
2652                 trD_s2(IJK,M) = ZERO ! initialize the totalizer
2653                 DO I1 = 1,3
2654                    DO I2 = 1,3
2655                       trD_s2(IJK,M) = trD_s2(IJK,M) + D_sM(I1,I2)*D_sM(I1,I2)
2656                    ENDDO
2657                 ENDDO
2658     
2659     ! use this fact to prevent underflow during theta calculation
2660                 IF (trD_s2(IJK,M) == zero) trD_s_C(IJK,M) = zero
2661     
2662     ! Frictional-flow stress tensor
2663     ! Needed by schaeffer or friction
2664     ! Calculate the second invariant of the deviator of D_s
2665                 I2_devD_s(IJK) = ( (D_sM(1,1)-D_sM(2,2))**2 + &
2666                                    (D_sM(2,2)-D_sM(3,3))**2 + &
2667                                    (D_sM(3,3)-D_sM(1,1))**2 )/6.&
2668                      + D_sM(1,2)**2 + D_sM(2,3)**2 + D_sM(3,1)**2
2669     
2670                 IF (V_EX /= ZERO) CALL CALC_BOYLE_MASSOUDI_STRESS(IJK,M,D_sM)
2671     
2672     ! Quantities for iddir-arastoopour theory: the trace of D_sm dot D_sl
2673     ! is required
2674     ! -------------------------------------------------------------------
2675                 IF (KT_TYPE_ENUM == IA_2005) THEN
2676                    DO L = 1,SMAX
2677                       IF (L .NE. M) THEN
2678                          IF (L > M) THEN  !done because trD_s2_ip(IJK,M,L) is symmetric, sof.
2679     
2680     ! Velocity derivatives (gradient and rate of strain tensor) for Mth
2681     ! solids phase at i, j, k
2682                             CALL CALC_DERIV_VEL_SOLIDS(IJK, L, DelV_sL, D_sL)
2683     
2684     ! Calculate trace of the D_sl dot D_sm
2685     ! (normal matrix multiplication)
2686                             trD_s2_ip(IJK,M,L) = ZERO ! initialize
2687                             DO I1 = 1,3
2688                                DO I2 = 1,3
2689                                   trD_s2_ip(IJK,M,L) = trD_s2_ip(IJK,M,L)+&
2690                                      D_sL(I1,I2)*D_sM(I1,I2)
2691                                ENDDO
2692                             ENDDO
2693     
2694                          ELSE   ! elseif (m<=m)
2695                             trD_s2_ip(IJK,M,L) = trD_s2_ip(IJK,L,M)
2696                          ENDIF ! end if/else (l>m)
2697                       ELSE   ! elseif (L=M)
2698                          trD_s2_ip(IJK,M,M) = trD_s2(IJK,M)
2699                       ENDIF   ! end if/else (m.NE.l)
2700                    ENDDO   ! end do (l=1,smax)
2701                 ENDIF   ! endif (kt_type = IA theory)
2702     
2703     
2704              ENDIF   ! end if (fluid_at)
2705           ENDDO   ! end outer IJK loop
2706     !!$omp end parallel do
2707     
2708           IF (SHEAR) THEN
2709              call remove_shear(M)
2710              IF (KT_TYPE_ENUM == IA_2005) THEN
2711                 DO L = 1,SMAX
2712                    IF (L .NE. M) THEN
2713                       CALL remove_shear(L)
2714                    ENDIF
2715                 ENDDO
2716              ENDIF
2717           ENDIF
2718     
2719           RETURN
2720           END SUBROUTINE INIT1_MU_S
2721     
2722     
2723     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2724     !                                                                      C
2725     !  Subroutine: calc_boyle_massoudi_stress                              C
2726     !  Purpose: Define quantities needed for viscosity calculations        C
2727     !  pertaining to the boyle-massoudi stress term.                       C
2728     !  Notes: this is applicable within the algebraic granular energy      C
2729     !  equation only.                                                      C
2730     !                                                                      C
2731     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2732           SUBROUTINE CALC_BOYLE_MASSOUDI_STRESS(IJK, M, D_S)
2733     
2734     ! Modules
2735     !---------------------------------------------------------------------//
2736           USE fldvar, only: ep_g, ep_s
2737           USE param1, only: zero, one
2738           USE visc_s, only: ep_star_array
2739           USE visc_s, only: trm_s, trdM_s
2740     
2741           USE geometry, only: odx_e, ody_n, odz_t, ox, x
2742           USE indices, only: i_of, j_of, k_of
2743           USE indices, only: im1, jm1, km1
2744           USE functions, only: west_of, east_of, south_of, north_of
2745           USE functions, only: top_of, bottom_of, fluid_at
2746           IMPLICIT NONE
2747     
2748     ! Local variables
2749     !---------------------------------------------------------------------//
2750     ! ijk index
2751           INTEGER, INTENT(IN) :: IJK
2752     ! solids phase index
2753           INTEGER, INTENT(IN) ::M
2754     ! rate of strain tensor at i, j, k
2755           DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: D_S
2756     
2757     ! Local variables
2758     !---------------------------------------------------------------------//
2759     ! Indices
2760           INTEGER :: I, J, K, IM, JM, KM
2761           INTEGER :: IJKW, IJKE, IJKS, IJKN, IJKB, IJKT
2762     ! Solids volume fraction gradient tensor
2763           DOUBLE PRECISION :: M_s(3,3)
2764     ! d(EP_sm)/dX
2765           DOUBLE PRECISION :: DEP_soDX
2766     ! d(EP_sm)/dY
2767           DOUBLE PRECISION :: DEP_soDY
2768     ! d(EP_sm)/XdZ
2769           DOUBLE PRECISION :: DEP_soXDZ
2770     ! Local DO-LOOP counters
2771           INTEGER :: I1, I2
2772     !---------------------------------------------------------------------//
2773     
2774           I = I_OF(IJK)
2775           J = J_OF(IJK)
2776           K = K_OF(IJK)
2777           IM = Im1(I)
2778           JM = Jm1(J)
2779           KM = Km1(K)
2780           IJKW  = WEST_OF(IJK)
2781           IJKE  = EAST_OF(IJK)
2782           IJKS  = SOUTH_OF(IJK)
2783           IJKN  = NORTH_OF(IJK)
2784           IJKB  = BOTTOM_OF(IJK)
2785           IJKT  = TOP_OF(IJK)
2786     
2787           IF(EP_g(IJK) .GE. EP_star_array(IJK)) THEN
2788              DEP_soDX  = ( EP_s(IJKE, M) - EP_s(IJK, M) ) * oDX_E(I)&
2789                  * ( ONE / ( (oDX_E(IM)/oDX_E(I)) + ONE ) ) +&
2790                  ( EP_s(IJK, M) - EP_s(IJKW, M) ) * oDX_E(IM)&
2791                  * ( ONE / ( (oDX_E(I)/oDX_E(IM)) + ONE ) )
2792                 DEP_soDY  = ( EP_s(IJKN, M) - EP_s(IJK, M) ) * oDY_N(J)&
2793                  * ( ONE / ( (oDY_N(JM)/oDY_N(J)) + ONE ) ) +&
2794                  ( EP_s(IJK, M) - EP_s(IJKS, M) ) * oDY_N(JM)&
2795                  * ( ONE / ( (oDY_N(J)/oDY_N(JM)) + ONE ) )
2796              DEP_soXDZ  = (( EP_s(IJKT, M) - EP_s(IJK, M) )&
2797                  * oX(I)*oDZ_T(K)&
2798                  * ( ONE / ( (oDZ_T(KM)/oDZ_T(K)) + ONE ) ) +&
2799                  ( EP_s(IJK, M) - EP_s(IJKB, M) ) * oX(I)*oDZ_T(KM)&
2800                  * ( ONE / ( (oDZ_T(K)/oDZ_T(KM)) + ONE ) ) ) /&
2801                  X(I)
2802     
2803               M_s(1,1) = DEP_soDX * DEP_soDX
2804               M_s(1,2) = DEP_soDX * DEP_soDY
2805               M_s(1,3) = DEP_soDX * DEP_soXDZ
2806               M_s(2,1) = DEP_soDX * DEP_soDY
2807               M_s(2,2) = DEP_soDY * DEP_soDY
2808               M_s(2,3) = DEP_soDY * DEP_soXDZ
2809               M_s(3,1) = DEP_soDX * DEP_soXDZ
2810               M_s(3,2) = DEP_soDY * DEP_soXDZ
2811               M_s(3,3) = DEP_soXDZ * DEP_soXDZ
2812     
2813               trM_s(IJK) = M_s(1,1) + M_s(2,2) + M_s(3,3)
2814               trDM_s(IJK) = ZERO
2815               DO I1 = 1,3
2816                  DO I2 = 1,3
2817                     trDM_s(IJK) = trDM_s(IJK) + D_s(I1,I2)*M_s(I1,I2)
2818                  ENDDO
2819               ENDDO
2820           ENDIF   ! end if (ep_g >=ep_star_array)
2821     
2822           RETURN
2823           END SUBROUTINE CALC_BOYLE_MASSOUDI_STRESS
2824