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