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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module: kintheory                                                   C
4     !  Purpose: Common block containing constants, variables, functions    C
5     !  used by various kinetic theory models                               C
6     !                                                                      C
7     !  Literature/Document References:                                     C
8     !                                                                      C
9     !  Garzo, V., Tenneti, S., Subramaniam, S., and Hrenya, C. M.,         C
10     !      "Enskog kinetic theory for monodisperse gas-solid flows", JFM,  C
11     !      Vol. 712, 2012, pp. 129-168                                     C
12     !                                                                      C
13     !                                                                      C
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
15     
16           MODULE kintheory
17     
18     
19     ! coefficient terms needed for stress (in addition to
20     ! mu_s and lambda_s which are defined in visc_s_mod, allocated
21     ! in allocate_arrays and initialized in set_constprop; and
22     ! P_s which is defined in fldvar_mod, allocated in
23     ! allocate_arrays, and initialized in init_fvars)
24     !     stress term with gradient in particle M velocity
25           DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: MU_sM_ip
26     !     stress term with gradient in particle L velocity
27           DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: MU_sL_ip
28     !     stress term with trace in particle M velocity
29           DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: XI_sM_ip
30     !     stress term with trace in particle L velocity
31           DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: XI_sL_ip
32     
33     
34     ! coefficient terms needed for momentum source (in addition to
35     ! F_SS which is defined in drag_mod, allocated in allocate_arrays
36     ! and initialized in set_constprop)
37     !     momentum source term with gradient in number density
38           DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: Fnu_s_ip
39     !     momentum source term with gradient in mixture temperature or
40     !     with the gradient in temperature of species M
41           DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: FT_sM_ip
42     !     momentum source term with gradient in temperature of species L
43           DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: FT_sL_ip
44     
45     
46     ! coefficient terms needed for heat flux (in addition to
47     ! kth_s which is defined in set_constprop, allocated
48     ! allocate_arrays and initialized in init_fvars)
49     !     heat flux term with gradient in granular temperature of species L
50           DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: Kth_sL_ip
51     !     heat flux term with gradient in number density of species M
52           DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: Knu_sM_ip
53     !     heat flux term with gradient in number density of species L
54           DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: Knu_sL_ip
55     !     heat flux term with velocity difference
56           DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: Kvel_s_ip
57     
58     
59     ! coefficient terms needed for energy dissipation
60     !     energy dissipation with difference in species granular
61     !     temperature: transfer between solid solid phases
62           DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: ED_ss_ip
63     !     energy dissipation term
64           DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: EDT_s_ip
65     !     energy dissipation with divergence of velocity of species M
66           DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: EDvel_sM_ip
67     !     energy dissipation with divergence of velocity of species L
68           DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: EDvel_sL_ip
69     !     coefficient A2, xsi used in multiple places in GTSH theory
70           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: A2_gtsh
71           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xsi_gtsh
72     
73     ! Solids source terms needed for Iddir & Arastoopour (2005)
74     ! kinetic theory model
75           DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE ::  KTMOM_U_s
76           DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE ::  KTMOM_V_s
77           DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE ::  KTMOM_W_s
78     
79     ! parameter in the theory of GTSH that is related to length scale
80     ! of lubrication effects. For details see GTSH, 2012.
81           DOUBLE PRECISION, PARAMETER :: EpM = 0.01d0
82     
83     
84     
85           CONTAINS
86     
87     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
88     !                                                                      C
89     !  Purpose: Calculate the magnitude of the gas-solids relative         C
90     !  velocity at i, j, k                                                 C
91     !                                                                      C
92     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
93           DOUBLE PRECISION FUNCTION KT_RVEL(IJK, m)
94     
95     ! Modules
96     !---------------------------------------------------------------------//
97           use fldvar, only: u_s, v_s, w_s
98           use fldvar, only: u_g, v_g, w_g
99           use run, only: shear
100           use vshear, only: vsh
101           use indices, only: i_of
102           use functions, only: im_of, jm_of, km_of
103           use functions, only: fluid_at
104           use fun_avg, only: avg_x_e, avg_y_n, avg_z_t
105           IMPLICIT NONE
106     
107     ! Dummy arguments
108     !---------------------------------------------------------------------//
109     ! ijk index
110           INTEGER, INTENT(IN) :: ijk
111     ! solids phase index
112           INTEGER, INTENT(IN) :: M
113     
114     ! Local variables
115     !---------------------------------------------------------------------//
116     ! cell indices
117           INTEGER :: I, IMJK, IJMK, IJKM
118     ! Cell center value of solids and gas velocities
119           DOUBLE PRECISION :: USCM, VSCM, WSCM, &
120                               UGC, VGC, WGC
121     ! y-component of velocity with 'shear' applied
122           DOUBLE PRECISION :: vs_j, vs_jm
123     !---------------------------------------------------------------------//
124           I = I_OF(IJK)
125           IMJK  = IM_OF(IJK)
126           IJMK  = JM_OF(IJK)
127           IJKM  = KM_OF(IJK)
128     
129     ! Awkward here but captures what was done in early version of
130     ! calc_mu_s. Otherwise any call to this routine must be done with
131     ! 'shear' already applied to v_s.
132           vs_j = V_S(IJK,M)
133           vs_jm = V_S(IJMK,M)
134           IF (SHEAR) THEN
135              vs_j = V_S(IJK,M)+VSH(IJK)
136              IF(FLUID_AT(IJMK)) THEN
137                 vs_jm = V_S(IJMK,M)+VSH(IJMK)
138              ELSE
139                 vs_jm = v_s(IJMK,M)
140              ENDIF
141           ENDIF
142     
143     ! Start calculation for relative velocity
144     ! Calculate velocity components at i, j, k
145           UGC = AVG_X_E(U_G(IMJK),U_G(IJK),I)
146           VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
147           WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
148     
149           USCM = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
150           VSCM = AVG_Y_N(vs_jm, vs_j)
151           WSCM = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M))
152     
153     ! Magnitude of gas-solids relative velocity
154           KT_RVEL = SQRT((UGC - USCM)**2 + (VGC - VSCM)**2 + &
155                          (WGC - WSCM)**2)
156     
157           RETURN
158           END FUNCTION KT_RVEL
159     
160     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
161     !                                                                      C
162     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
163           DOUBLE PRECISION FUNCTION KT_COS_THETA(ijk, M)
164     
165     ! Modules
166     !---------------------------------------------------------------------//
167           USE param1, only: zero, small_number
168           USE fldvar, only: u_g, v_g, w_g
169           USE fldvar, only: u_s, v_s, w_s
170           USE fldvar, only: ep_s
171           USE run, only: shear
172           USE vshear, only: vsh
173           USE toleranc, only: zero_ep_s
174           USE indices, only: i_of
175           USE functions, only: im_of, jm_of, km_of
176           USE functions, only: fluid_at
177           USE fun_avg, only: avg_x_e, avg_y_n, avg_z_t
178           IMPLICIT NONE
179     
180     ! Dummy arguments
181     !---------------------------------------------------------------------//
182     ! ijk index
183           INTEGER, INTENT(IN) :: ijk
184     ! solids phase index
185           INTEGER, INTENT(IN) :: M
186     
187     ! Local variables
188     !---------------------------------------------------------------------//
189     ! cell indices
190           INTEGER :: I, IMJK, IJMK, IJKM
191     ! Cell center value of solids and gas velocities
192           DOUBLE PRECISION :: USCM, VSCM, WSCM, &
193                               UGC, VGC, WGC
194           DOUBLE PRECISION :: vs_j, vs_jm, speed, rvel
195     !---------------------------------------------------------------------//
196     
197           I = I_OF(IJK)
198           IMJK  = IM_OF(IJK)
199           IJMK  = JM_OF(IJK)
200           IJKM  = KM_OF(IJK)
201     
202     ! Awkward here but captures what was done in early version of
203     ! calc_mu_s. Otherwise any call to this routine must be done with
204     ! 'shear' already applied to v_s.
205           vs_j = V_S(IJK,M)
206           vs_jm = V_S(IJMK,M)
207           IF (SHEAR) THEN
208              vs_j = V_S(IJK,M)+VSH(IJK)
209              IF(FLUID_AT(IJMK)) THEN
210                 vs_jm = V_S(IJMK,M)+VSH(IJMK)
211              ELSE
212                 vs_jm = v_s(IJMK,M)
213              ENDIF
214           ENDIF
215     
216           UGC = AVG_X_E(U_G(IMJK),U_G(IJK),I)
217           VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
218           WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
219     
220           USCM = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
221           VSCM = AVG_Y_N(vs_jm, vs_j)
222           WSCM = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M))
223     
224           RVEL =  SQRT((UGC - USCM)**2 + (VGC - VSCM)**2 + &
225                        (WGC - WSCM)**2)
226     
227           SPEED = SQRT(USCM**2+VSCM**2+WSCM**2)
228     
229     ! motion viewed by the particles (crossing trajectory effect)
230           IF(SPEED > Small_Number .AND. RVEL > Small_Number .AND. &
231              EP_S(IJK,M) > ZERO_EP_S) THEN
232              KT_Cos_Theta = ( (UGC-USCM)*USCM + (VGC-VSCM)*VSCM + &
233                            (WGC-WSCM)*WSCM )/(RVEL*SPEED)
234           ELSE
235              KT_Cos_Theta = ZERO
236           ENDIF
237     
238           RETURN
239           END FUNCTION  KT_COS_THETA
240     
241     
242     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
243     !                                                                      C
244     !  Purpose: Calculate a single particle drag coefficient/ep_s, which   C
245     !  is used in evaluting limiting values of certain granular kinetic    C
246     !  theory terms                                                        C
247     !                                                                      C
248     !  Comments: This is currently based on wen-yu but in the the future   C
249     !  we may want to make this consistent with drag_gs (that is, replace  C
250     !  this function call and assign appropriate variable within drag_gs   C
251     !                                                                      C
252     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
253           DOUBLE PRECISION FUNCTION KT_DGA(IJK, m)
254     
255     ! Modules
256     !---------------------------------------------------------------------//
257           use param1, only: zero, one, small_number, large_number
258           use fldvar, only: rop_g
259           use fldvar, only: d_p
260           use physprop, only: mu_g
261           IMPLICIT NONE
262     
263     ! Dummy arguments
264     !---------------------------------------------------------------------//
265     ! ijk index
266           INTEGER, INTENT(IN) :: IJK
267     ! solids phase index
268           INTEGER, INTENT(IN) :: M
269     
270     ! Local variables
271     !---------------------------------------------------------------------//
272     ! single particle drag coefficient, reynolds number
273           DOUBLE PRECISION :: C_d, Re
274     ! local value for relative velocity
275           DOUBLE PRECISION :: rvel
276     !---------------------------------------------------------------------//
277     
278     ! initialization
279           kt_dga = zero
280           rvel = zero
281           RVEL = KT_RVEL(IJK, M)
282     
283     ! Defining single particle drag coefficient dgA based on Wen-Yu
284     ! correlation
285           RE = D_p(IJK,M)*RVEL*ROP_G(IJK)/(MU_G(IJK) + SMALL_NUMBER)
286           IF(RE .LE. 1000.d0)THEN
287              C_d = (24.d0/(Re+SMALL_NUMBER)) * &
288                 (ONE + 0.15d0 * Re**0.687D0)
289           ELSE
290              C_d = 0.44d0
291           ENDIF
292     
293     ! dga_s is local to this routine as is lrvel
294           kt_dga = 0.75d0 * C_d * RVEL * ROP_g(IJK) / D_p(IJK,M)
295     
296     ! set value for 1st iteration and 1st time step
297           IF(RVEL == ZERO) kt_dga = LARGE_NUMBER
298     
299           RETURN
300           END FUNCTION KT_DGA
301     
302     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
303     !                                                                      C
304     !  Subroutine: CALC_IA_ENERGY_DISSIPATION_SS                           C
305     !                                                                      C
306     !  Purpose: Implement kinetic theory of Iddir & Arastoopour (2005)     C
307     !     for calculation of source terms in granular energy equation      C
308     !                                                                      C
309     !  Author: Janine E. Galvin, Univeristy of Colorado                    C
310     !                                                                      C
311     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
312           SUBROUTINE CALC_IA_ENERGY_DISSIPATION_SS(M)
313     
314     ! Modules
315     !---------------------------------------------------------------------//
316           USE constant, only: pi
317           USE constant, only: C_e
318           USE fldvar, only: ro_s, rop_s, d_p
319           USE fldvar, only: theta_m
320           USE physprop, only: mmax
321           USE rdf, only: g_0
322           USE param1, only: zero
323     
324           USE functions, only: fluid_at, funlm
325           USE compar, only: ijkstart3, ijkend3
326           IMPLICIT NONE
327     
328     ! Dummy arguments
329     !---------------------------------------------------------------------//
330     ! Solids phase index
331           INTEGER, INTENT(IN) :: M
332     
333     ! Local variables
334     !---------------------------------------------------------------------//
335     ! Index
336           INTEGER :: IJK
337     ! Solids phase index
338           INTEGER :: L
339     ! Index for storing solids-solids drag coefficients
340     ! in the upper triangle of the matrix
341           INTEGER :: LM
342     ! variables for IA theory
343           DOUBLE PRECISION :: ED_common_term
344           DOUBLE PRECISION :: EDvel_sL, EDvel_sM
345           DOUBLE PRECISION :: M_PM, M_PL, MPSUM, NU_PL, NU_PM, D_PM, &
346                               D_PL, DPSUMo2
347           DOUBLE PRECISION :: Ap_lm, Dp_lm, R1p_lm, R10p_lm, R3p_lm, &
348                               R4p_lm, R5p_lm, Bp_lm
349     !---------------------------------------------------------------------//
350     
351           DO IJK = ijkstart3, ijkend3
352               IF ( FLUID_AT(IJK) ) THEN
353     
354                  D_PM = D_P(IJK,M)
355                  M_PM = (PI/6.d0)*(D_PM**3)*RO_S(IJK,M)
356                  NU_PM = ROP_S(IJK,M)/M_PM
357     
358                  DO L = 1, MMAX
359                     LM = FUNLM(L,M)
360                     D_PL = D_P(IJK,L)
361                     M_PL = (PI/6.d0)*(D_PL**3)*RO_S(IJK,L)
362     
363                     MPSUM = M_PM + M_PL
364                     DPSUMo2 = (D_PM+D_PL)/2.d0
365                     NU_PL = ROP_S(IJK,L)/M_PL
366     
367                     ED_common_term = (3.d0/4.d0)*(DPSUMo2*DPSUMo2)*&
368                        (1.d0+C_E)*G_0(IJK,M,L)*NU_PM*NU_PL*(M_PM*&
369                        M_PL/MPSUM)*((M_PM*M_PL)**1.5)
370     
371                     IF (M .eq. L) THEN
372                        Ap_lm = MPSUM/(2.d0)
373                        Dp_lm = M_PL*M_PM/(2.d0*MPSUM)
374                        R1p_lm = 1.d0/( (Ap_lm**1.5)*(Dp_lm**3) )
375                        R3p_lm = 1.d0/( (Ap_lm**1.5)*(Dp_lm**3.5) )
376     
377     ! Dissipation associated with the difference in temperature
378     ! (e.g. interphase transfer term).  For the case case (L=M)
379     ! the term ED_s * (TL-TM) cancels. Therefore, explicity set
380     ! the term to zero for this case.
381                        ED_ss_ip(IJK,LM) = ZERO
382     
383     ! Dissipation associated with temperature
384                        EDT_s_ip(IJK,M,L) = -ED_common_term* (1.d0-C_E)*&
385                           (M_PL/MPSUM)*(DSQRT(PI)/6.d0)*R1p_lm*&
386                           (Theta_m(IJK,M)**1.5)
387     
388     ! Dissipation associated with divergence of
389     ! velocity of solid phase L: do not explicity include
390     ! terms that cancel when EDvel_sL is summed with EDvel_sM
391                        EDvel_sL = ED_common_term*((1.d0-C_E)*(M_PL/MPSUM)*&
392                           (DPSUMo2*PI/48.d0)*(M_PM*M_PL/MPSUM)*R3p_lm)
393                           EDvel_sL_ip(IJK,M,L) = EDvel_sL*(Theta_m(IJK,L))
394     
395     ! Dissipation associated with divergence of
396     ! velocity of solid phase M: do not explicity include
397     ! terms that cancel when EDvel_sL is summed with EDvel_sM
398     ! commented by sof, no need to re-do computation
399                        !EDvel_sM = ED_common_term*((1.d0-C_E)*(M_PL/MPSUM)*&
400                        !   (DPSUMo2*PI/48.d0)*(M_PM*M_PL/MPSUM)*R3p_lm)
401                        !EDvel_sM_ip(IJK,M,L) = EDvel_sM*(Theta_m(IJK,M)) ! M=L
402     
403                        EDvel_sM_ip(IJK,M,L) =  EDvel_sL_ip(IJK,M,L)
404     
405                     ELSE
406     
407                        Ap_lm = (M_PM*Theta_m(IJK,L)+M_PL*Theta_m(IJK,M))/&
408                           2.d0
409                        Bp_lm = (M_PM*M_PL*(Theta_m(IJK,L)-&
410                           Theta_m(IJK,M) ))/(2.d0*MPSUM)
411                        Dp_lm = (M_PL*M_PM*(M_PM*Theta_m(IJK,M)+M_PL*&
412                           Theta_m(IJK,L) ))/(2.d0*MPSUM*MPSUM)
413     
414                        R1p_lm = (1.d0/((Ap_lm**1.5)*(Dp_lm**3)))+ &
415                           ((9.d0*Bp_lm*Bp_lm)/(Ap_lm**2.5 * Dp_lm**4))+&
416                           ((30.d0*Bp_lm**4)/(2.d0*Ap_lm**3.5 * Dp_lm**5))
417     
418                        R3p_lm = (1.d0/((Ap_lm**1.5)*(Dp_lm**3.5)))+&
419                           ((21.d0*Bp_lm*Bp_lm)/(2.d0 * Ap_lm**2.5 * Dp_lm**4.5))+&
420                           ((315.d0*Bp_lm**4)/(8.d0 * Ap_lm**3.5 *Dp_lm**5.5))
421     
422                        R4p_lm = (3.d0/( Ap_lm**2.5 * Dp_lm**3.5))+&
423                           ((35.d0*Bp_lm*Bp_lm)/(2.d0 * Ap_lm**3.5 * Dp_lm**4.5))+&
424                           ((441.d0*Bp_lm**4)/(8.d0 * Ap_lm**4.5 * Dp_lm**5.5))
425     
426                        R5p_lm = (1.d0/(Ap_lm**2.5 * Dp_lm**3))+ &
427                           ((5.d0*Bp_lm*Bp_lm)/(Ap_lm**3.5 * Dp_lm**4))+&
428                           ((14.d0*Bp_lm**4)/( Ap_lm**4.5 * Dp_lm**5))
429     
430                        R10p_lm = (1.d0/(Ap_lm**2.5 * Dp_lm**2.5))+&
431                           ((25.d0*Bp_lm*Bp_lm)/(2.d0* Ap_lm**3.5 * Dp_lm**3.5))+&
432                           ((1225.d0*Bp_lm**4)/(24.d0* Ap_lm**4.5 * Dp_lm**4.5))
433     
434     ! Dissipation associated with the difference in temperature
435     ! (e.g. interphase transfer term). to solved using PEA
436                        ED_ss_ip(IJK,LM) = ED_common_term*DSQRT(PI)*(M_PM*M_PL/&
437                           (2.d0*MPSUM))*R5p_lm*( (Theta_M(IJK,M)*Theta_M(IJK,L))**3 )
438     
439     ! Dissipation associated with temperature
440                        EDT_s_ip(IJK,M,L) = -ED_common_term* (1.d0-C_E)*(M_PL/MPSUM)*&
441                           (DSQRT(PI)/6.d0)*R1p_lm*( (Theta_m(IJK,M)*Theta_m(IJK,L))**3 )
442     
443     ! Dissipation associated with divergence of
444     ! velocity of solid phase L
445                        EDvel_sL = ED_common_term*( ((3.d0*DPSUMo2*PI/40.d0)*M_PL*&
446                           R10p_lm)+( (DPSUMo2*PI/4.d0)*(M_PM*M_PL/MPSUM)*Bp_lm*&
447                           R4p_lm)-( (1.d0-C_E)*(M_PL/MPSUM)*(DPSUMo2*PI/16.d0)*&
448                           M_PL*Bp_lm*R4p_lm)+( (1.d0-C_E)*(M_PL/MPSUM)*&
449                           (DPSUMo2*PI/48.d0)*(M_PM*M_PL/MPSUM)*R3p_lm))
450     
451                        EDvel_sL_ip(IJK,M,L) = EDvel_sL*( Theta_m(IJK,M)**3.5 *&
452                           Theta_m(IJK,L)**2.5 )
453     
454     ! Dissipation associated with divergence of
455     ! velocity of solid phase M
456                        EDvel_sM = ED_common_term*( (-(3.d0*DPSUMo2*PI/40.d0)*M_PM*&
457                           R10p_lm)+( (DPSUMo2*PI/4.d0)*(M_PM*M_PL/MPSUM)*Bp_lm*&
458                           R4p_lm)+( (1.d0-C_E)*(M_PL/MPSUM)*(DPSUMo2*PI/16.d0)*&
459                            M_PM*Bp_lm*R4p_lm)+( (1.d0-C_E)*(M_PL/MPSUM)*&
460                           (DPSUMo2*PI/48.d0)*(M_PM*M_PL/MPSUM)*R3p_lm))
461     
462                        EDvel_sM_ip(IJK,M,L) = EDvel_sM*( Theta_m(IJK,M)**2.5 *&
463                           Theta_m(IJK,L)**3.5 )
464     
465                     ENDIF   ! end if/else (m.eq.l)
466     
467                  ENDDO   ! end do (l=1,mmax)
468     
469              ENDIF   ! end if(fluid_at)
470           ENDDO   ! end do ijk
471     
472           RETURN
473           END SUBROUTINE CALC_IA_ENERGY_DISSIPATION_SS
474     
475     
476     
477     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
478     !                                                                      C
479     !  Subroutine: CALC_GD_99_ENERGY_DISSIPATION_SS                        C
480     !                                                                      C
481     !  Purpose: Implement kinetic theory of Garzo & Dufty (1999)           C
482     !  for calculation of source terms in granular energy equation         C
483     !                                                                      C
484     !  Author: Janine E. Galvin                                            C
485     !                                                                      C
486     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
487           SUBROUTINE CALC_GD_99_ENERGY_DISSIPATION_SS(M)
488     
489     ! Modules
490     !---------------------------------------------------------------------//
491           USE constant, only: pi
492           USE constant, only: C_e
493           USE fldvar, only: rop_s, d_p
494           USE fldvar, only: theta_m, ep_s
495           USE rdf, only: g_0
496     
497           USE functions, only: fluid_at
498           USE compar, only: ijkstart3, ijkend3
499           IMPLICIT NONE
500     
501     ! Dummy arguments
502     !---------------------------------------------------------------------//
503     ! Solids phase index
504           INTEGER, INTENT(IN) :: M
505     
506     ! Local variables
507     !---------------------------------------------------------------------//
508     ! Indices
509           INTEGER :: IJK
510     ! variables for GD model
511           DOUBLE PRECISION :: press_star, c_star, zeta0_star, &
512                               nu_gamma_star, &
513                               lambda_num, cd_num, zeta1
514           DOUBLE PRECISION :: D_PM, EP_SM
515           DOUBLE PRECISION :: nu0, Chi
516     !---------------------------------------------------------------------//
517     
518           DO IJK = ijkstart3, ijkend3
519               IF ( FLUID_AT(IJK) ) THEN
520     
521     ! Note: k_boltz = M_PM
522     
523     ! local aliases
524                  Chi = G_0(IJK,M,M)
525                  EP_SM = EP_s(IJK,M)
526                  D_PM = D_P(IJK,M)
527     
528     !            nu0=p_k/eta0
529                  nu0 = (96.d0/5.d0)*(EP_SM/D_PM)*DSQRT(Theta_m(IJK,M)/PI)
530     
531                  press_star = 1.d0 + 2.d0*(1.d0+C_E)*EP_SM*Chi
532     
533                  c_star = 32.0d0*(1.0d0 - C_E)*(1.d0 - 2.0d0*C_E*C_E) &
534                     / (81.d0 - 17.d0*C_E + 30.d0*C_E*C_E*(1.0d0-C_E))
535     
536                  zeta0_star = (5.d0/12.d0)*Chi*(1.d0 - C_E*C_E) &
537                     * (1.d0 + (3.d0/32.d0)*c_star)
538     
539                  nu_gamma_star = ((1.d0+C_E)/48.d0)*Chi*(128.d0-96.d0*C_E + &
540                     15.d0*C_E*C_E-15.d0*C_E*C_E*C_E+ (c_star/64.d0) * (15.d0* &
541                     C_E*C_E*C_E-15.d0*C_E*C_E+498.d0*C_E-434.d0))
542     
543                  lambda_num = (3.d0/8.d0)*( (1.d0-C_E)*(5.d0*C_E*C_E+4.d0*C_E-1.d0)+ &
544                     (c_star/12.d0)*(159.d0*C_E+3.d0*C_E*C_E-19.d0*C_E- &
545                     15.d0*C_E*C_E*C_E) ) * (1.d0+C_E)
546     
547     ! does not include factor of 1/nu0.  the factor 1/nu0 in cD will cancel
548     ! with the factor of nu0 used to obtain zeta1 from zeta1_star
549                  cd_num = ( (4.d0/15.d0)*lambda_num*EP_SM*Chi + &
550                     (press_star-1.d0)*(2.d0/3.d0-C_E)*c_star ) / &
551                     ( 0.5d0*zeta0_star+nu_gamma_star + (5.d0*c_star/64.d0) * &
552                     (1.d0+(3.d0*c_star/64.d0))*Chi*(1.d0-C_E*C_E))
553     
554     ! does not include factor of 1/nu0 in the first term.  the factor 1/nu0
555     ! will cancel with the factor of nu0 used to obtain zeta1 from zeta1_star
556     ! zeta1 = nu0*zeta1_star
557                  zeta1 = -(1.d0-C_E)*(press_star-1.d0) + (5.d0/32.d0) * &
558                     (1.d0-C_E*C_E)*(1.d0+(3.d0*c_star/64.d0))*Chi*cd_num
559     
560     ! in the energy equation this term is multiplied by 3/2*n*kboltz*T
561     ! leave multiplication of theta for source routine
562                  EDvel_sM_ip(IJK,M,M) = (3.d0/2.d0)*ROP_s(IJK,M)*zeta1
563     
564     ! in the energy equation this term is multiplied by 3/2*n*kboltz*T
565     ! leave multiplication of theta for source routine
566                  EDT_s_ip(IJK,M,M) = (3.d0/2.d0)*ROP_s(IJK,M)*nu0*zeta0_star
567     
568              ENDIF   ! end if (fluid_at)
569           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
570     
571           RETURN
572           END SUBROUTINE CALC_GD_99_ENERGY_DISSIPATION_SS
573     
574     
575     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
576     !                                                                      C
577     !  Subroutine: CALC_GTSH_ENERGY_DISSIPATION_SS                         C
578     !                                                                      C
579     !  Purpose: Implement kinetic theory of Garzo, Tenneti, Subramaniam    C
580     !  Hrenya (2012) for calculation of source terms in granular           C
581     !  energy equation                                                     C
582     !                                                                      C
583     !  Author: Sofiane Benyahia                                            C
584     !                                                                      C
585     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
586           SUBROUTINE CALC_GTSH_ENERGY_DISSIPATION_SS(M)
587     
588     ! Modules
589     !---------------------------------------------------------------------//
590           USE constant, only: pi
591           USE constant, only: C_e
592           USE fldvar, only: ep_s
593           USE fldvar, only: ro_g, rop_g
594           USE fldvar, only: ro_s, d_p
595           USE fldvar, only: theta_m
596           USE physprop, only: mu_g
597           USE rdf, only: g_0
598           USE param1, only: zero, one, small_number
599     
600           USE functions, only: fluid_at
601           USE compar, only: ijkstart3, ijkend3
602           IMPLICIT NONE
603     
604     ! Dummy arguments
605     !---------------------------------------------------------------------//
606     ! Solids phase index
607           INTEGER, INTENT(IN) :: M
608     
609     ! Local variables
610     !---------------------------------------------------------------------//
611     ! Indices
612           INTEGER :: IJK
613     ! variables
614           DOUBLE PRECISION :: D_PM, EP_SM, V_p, N_p, M_p
615           DOUBLE PRECISION :: nu0, Chi
616           DOUBLE PRECISION :: VREL
617           DOUBLE PRECISION :: Re_m, Re_T
618           DOUBLE PRECISION :: zeta_star, mu2_0, mu4_0, mu4_1
619           DOUBLE PRECISION :: omega, nu_j, rho_10, rho_11
620     
621     !---------------------------------------------------------------------//
622     
623           DO IJK = ijkstart3, ijkend3
624              IF ( FLUID_AT(IJK) ) THEN
625     
626     ! Local aliases
627                 Chi = G_0(IJK,M,M)
628                 EP_SM = EP_s(IJK,M)
629                 D_PM = D_P(IJK,M)
630                 V_p = pi*D_PM**3/6d0
631                 n_p = EP_SM/V_p
632                 M_p = V_p * ro_s(ijk,m)
633     
634                 nu0 = (96.d0/5.d0)*(EP_SM/D_PM)*DSQRT(Theta_m(IJK,M)/PI)
635     
636     ! First calculate the Re number (Re_m, Re_T) in eq. 3.1 in GTSH theory
637                 VREL = KT_RVEL(IJK, M)
638     ! Note: rop_g = ro_g * (1-phi)
639                 Re_m = D_PM*VREL*ROP_g(ijk)/Mu_g(ijk)
640                 Re_T = ro_g(ijk)*D_PM*dsqrt(theta_m(ijk,m)) / mu_g(ijk)
641     
642     ! Now calculate and store xsi (used in many spots) in eq. 8.2 in
643     ! GTSH theory.
644                 xsi_gtsh(ijk) = one/6d0*D_PM*VREL**2*(3d0*pi*mu_g(ijk)*&
645                    D_PM/M_p)**2 / dsqrt(pi*theta_m(ijk,m)) * &
646                    S_star(EP_SM, Chi)
647     
648     ! eq. (6.22) GTSH theory
649                 mu2_0 = dsqrt(2d0*pi) * Chi * (one-C_E**2)
650     
651     ! eq. (6.23) GTSH theory
652                 mu4_0 = (4.5d0+C_E**2) * mu2_0
653     
654     ! eq. (6.24) GTSH theory
655     !            mu4_1 = (6.46875d0+0.3125d0*C_E**2 + 2d0/(one-C_E)) * mu2_0
656     ! this is done to avoid /0 in case c_e = 1.0
657                 mu4_1 = (6.46875d0+0.9375d0*C_E**2)*mu2_0 + &
658                    2d0*dsqrt(2d0*pi)* Chi*(one+C_E)
659     
660                 A2_gtsh(ijk) = zero ! for EP_SM = zero
661                 if(EP_SM> small_number) then ! avoid singularity
662     ! Now calculate zeta_star in eq. 8.10 in GTSH theory
663                    zeta_star = 4.5d0*dsqrt(2d0*Pi)*&
664                       (ro_g(ijk)/ro_s(ijk,m))**2*Re_m**2 * &
665                       S_star(EP_SM,Chi) / (EP_SM*(one-EP_SM)**2 * Re_T**4)
666     
667     ! Now calculate important parameter A2. This is used in many transport
668     ! coefficients so that it's best to store it in an array instead of
669     ! having many function calls.
670                    A2_gtsh(ijk) = (5d0*mu2_0 - mu4_0) / (mu4_1 - 5d0* &
671                       (19d0/16d0*mu2_0 - 1.5d0*zeta_star))
672                 endif
673     
674     ! Now calculate first term rho_0 in cooling rate rho, eq (6.26) in GTSH
675     ! theory.
676     ! note that theta_(ijk,m) does't contain mass m, so that T/m = theta_m.
677     ! note that edt_s_ip will need to be multiplied by 3/2 rop_s(ijk,m) in
678     ! source_granular_energy to have the same meaning as in GD_99 theory.
679                 EDT_s_ip(ijk,M,M) = 4d0/3d0*dsqrt(pi)*(one-C_E**2)*Chi* &
680                    (one+0.1875d0*A2_gtsh(ijk))*n_p*D_PM**2*dsqrt(theta_m(ijk,m))
681     
682     ! Calculate the second term in cooling rate, eq (7.25) in GTSH theory.
683     ! First calculate eq (7.28, 7.29) of GTSH theory.
684                 omega = (one+C_E)*nu0*((one-C_E**2)*(5d0*C_E-one) - &
685                    A2_gtsh(ijk)/ 6d0 * (15d0*C_E**3-3d0*C_E**2+81d0*C_E-61d0))
686                 nu_j = (one+C_E)/192d0*Chi*nu0* &
687                    (241d0-177d0*C_E+30d0*C_E**2-30d0*C_E**3)
688     
689     ! Now calculate eq (7.26, 7.27) of GTSH theory. ! corrected by W. Fullmer
690                 rho_10 = 2d0*Chi*EP_SM*(C_E**2-one)
691                 rho_11 = 25d0/1024d0*EP_SM*Chi**2*(one-C_E**2)* &
692                    (one+3d0/128d0*A2_gtsh(ijk)) * (omega/10d0 - &
693                    (one+C_E)*nu0*(one/3d0-C_E)*A2_gtsh(ijk)/2d0) / &
694                    (nu_j + G_gtsh(EP_SM, Chi, IJK, M)/m_p + 1.5d0* &
695                    xsi_gtsh(ijk)/theta_m(ijk,m) -1.5d0*EDT_s_ip(ijk,M,M))
696     
697     ! note that EDvel_sM_ip will later need to be multiplied by 3/2 rop_s(ijk,m).
698                 EDvel_sM_ip(IJK,M,M) = rho_10 + rho_11
699     
700              ENDIF   ! end if (fluid_at)
701           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
702     
703           RETURN
704           END SUBROUTINE CALC_GTSH_ENERGY_DISSIPATION_SS
705     
706     
707     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
708     !                                                                      C
709     !  Functions: G_tsh, K_phi, R_d, S_star                                C
710     !                                                                      C
711     !  Purpose: Implement kinetic theory of Garzo, Tenneti, Subramaniam    C
712     !  Hrenya (2012) for calculation of source terms in granular           C
713     !  energy equation                                                     C
714     !                                                                      C
715     !  Author: Sofiane Benyahia                                            C
716     !                                                                      C
717     !  Comments:                                                           C
718     !     these functions are needed only for gtsh theory                  C
719     !     Function gamma, eq. (8.1) in GTSH theory                         C
720     !                                                                      C
721     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
722           DOUBLE PRECISION FUNCTION G_gtsh (EP_SM, Chi, IJK, M)
723     
724     ! Modules
725     !---------------------------------------------------------------------//
726           USE constant, only: pi
727           USE param1, only: one
728           USE physprop, only: mu_g
729           USE fldvar, only: ro_g
730           USE fldvar, only: d_p, theta_m
731           IMPLICIT NONE
732     
733     ! Dummy arguments
734     !---------------------------------------------------------------------//
735     ! solids volume fraction of phase m at ijk
736           DOUBLE PRECISION, INTENT(IN) :: EP_SM
737     ! radial distribution function of phase M at ijk
738           DOUBLE PRECISION, INTENT(IN) :: Chi
739     ! Index
740           INTEGER, INTENT(IN) :: IJK
741     ! Solids phase index. Note M should be equal to 1 since theory
742     ! valid for only mmax = 1.
743           INTEGER, INTENT(IN) :: M
744     
745     ! Local variables
746     !---------------------------------------------------------------------//
747           DOUBLE PRECISION :: Re_T
748           DOUBLE PRECISION :: Rdiss, RdissP
749     
750     !---------------------------------------------------------------------//
751     
752           if(EP_SM <= 0.1d0) then
753              RdissP = one+3d0*dsqrt(EP_SM/2d0)
754           else
755              RdissP = &
756              one + 3d0*dsqrt(EP_SM/2d0) + 135d0/64d0*EP_SM*dlog(EP_SM) + &
757              11.26d0*EP_SM*(one-5.1d0*EP_SM+16.57d0*EP_SM**2-21.77d0*    &
758              EP_SM**3) - EP_SM*Chi*dlog(epM)
759           endif
760     
761           Re_T = ro_g(ijk)*d_p(ijk,m)*dsqrt(theta_m(ijk,m)) / mu_g(ijk)
762     
763           Rdiss = RdissP + Re_T * K_phi(EP_SM)
764     
765           G_gtsh = 3d0*Pi*mu_g(ijk)*d_p(ijk,m)*Rdiss  ! eq. (8.1) in GTSH theory
766     
767           RETURN
768           END FUNCTION G_gtsh
769     
770     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
771     !                                                                      C
772     !  Comments:                                                           C
773     !     Function K(phi), eq. (2.7) in Yi/Zenk/Mitrano/Hrenya JFM (2013)  C
774     !                                                                      C
775     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
776           DOUBLE PRECISION FUNCTION K_phi (phi)
777           IMPLICIT NONE
778     
779     ! Dummy Arguments
780     !---------------------------------------------------------------------//
781     ! solids volume fraction
782           DOUBLE PRECISION, INTENT(IN) :: phi
783     !---------------------------------------------------------------------//
784     
785           K_phi = (0.096d0 + 0.142d0*phi**0.212d0) / (1d0-phi)**4.454d0
786           K_phi = 0.0d0 ! set to zero for compatibility with GTSH JFM (2012)
787           RETURN
788           END FUNCTION K_phi
789     
790     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
791     !                                                                      C
792     !  Comments:                                                           C
793     !     Function R_d, eq. (8.6) in GTSH theory                           C
794     !                                                                      C
795     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
796           DOUBLE PRECISION FUNCTION R_d (phi)
797           IMPLICIT NONE
798     
799     ! Dummy arguments
800     !---------------------------------------------------------------------//
801     ! solids volume fraction
802           DOUBLE PRECISION, INTENT(IN) :: phi
803     !---------------------------------------------------------------------//
804     
805           R_d = 1.0d0  ! this avoids singularity at phi = 0.0
806           if((phi > 1d-15) .and. (phi <= 0.4d0)) then
807             R_d = (1d0+3d0*dsqrt(phi/2d0)+135d0/64d0*phi*dlog(phi)+17.14d0*phi) / &
808                   (1d0+0.681d0*phi-8.48*phi**2+8.16d0*phi**3)
809           elseif(phi > 0.4d0) then
810             R_d = 10d0*phi/(1d0-phi)**3 + 0.7d0
811           endif
812           RETURN
813           END FUNCTION R_d
814     
815     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
816     !                                                                      C
817     !  Comments:                                                           C
818     !     Function S_Star(phi), eq. (8.3, 8.5) in GTSH theory              C
819     !                                                                      C
820     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
821           DOUBLE PRECISION FUNCTION S_star (phi, Chi)
822           IMPLICIT NONE
823     
824     ! Dummy arguments
825     !---------------------------------------------------------------------//
826     ! solids volume fraction
827           DOUBLE PRECISION, INTENT(IN) :: phi
828     ! radial distribution function
829           DOUBLE PRECISION, INTENT(IN) :: Chi
830     
831     !---------------------------------------------------------------------//
832     
833           S_star = 1.0d0
834           if(phi >= 0.1d0) &
835             S_star = R_d(phi)**2/(Chi*(1d0+3.5d0*dsqrt(phi)+5.9*phi))
836           RETURN
837           END FUNCTION S_Star
838     
839           END MODULE kintheory
840