File: /nfs/home/0/users/jenkins/mfix.git/model/calc_mu_s.f

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