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