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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_IA_ENERGY_DISSIPATION_SS                     C
4     !                                                                      C
5     !  Purpose: Implement kinetic theory of Iddir & Arastoopour (2005)     C
6     !     for calculation of source terms in granular energy equation      C
7     !                                                                      C
8     !  Author: Janine E. Galvin, Univeristy of Colorado                    C
9     !                                                                      C
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
11     
12           SUBROUTINE CALC_IA_ENERGY_DISSIPATION_SS(M)
13     
14     !-----------------------------------------------
15     ! Modules
16     !-----------------------------------------------
17           USE compar
18           USE constant
19           USE fldvar
20           USE functions
21           USE geometry
22           USE indices
23           USE kintheory
24           USE param
25           USE param1
26           USE physprop
27           USE rdf
28           USE run
29           USE toleranc
30           IMPLICIT NONE
31     !-----------------------------------------------
32     ! Dummy arguments
33     !-----------------------------------------------
34     ! Solids phase index
35           INTEGER, INTENT(IN) :: M
36     !-----------------------------------------------
37     ! Local variables
38     !-----------------------------------------------
39     ! Index
40           INTEGER :: IJK, I, J, K
41     
42     ! Solids phase index
43           INTEGER :: L
44     
45     ! Index for storing solids-solids drag coefficients
46     ! in the upper triangle of the matrix
47           INTEGER :: LM
48     
49     ! variables for IA theory
50           DOUBLE PRECISION :: ED_common_term
51           DOUBLE PRECISION :: EDvel_sL, EDvel_sM
52           DOUBLE PRECISION :: M_PM, M_PL, MPSUM, NU_PL, NU_PM, D_PM, &
53                               D_PL, DPSUMo2
54           DOUBLE PRECISION :: Ap_lm, Dp_lm, R1p_lm, R10p_lm, R3p_lm, &
55                               R4p_lm, R5p_lm, Bp_lm
56     !-----------------------------------------------
57     
58           DO IJK = ijkstart3, ijkend3
59               I = I_OF(IJK)
60               J = J_OF(IJK)
61               K = K_OF(IJK)
62     
63               IF ( FLUID_AT(IJK) ) THEN
64     
65                  D_PM = D_P(IJK,M)
66                  M_PM = (PI/6.d0)*(D_PM**3)*RO_S(IJK,M)
67                  NU_PM = ROP_S(IJK,M)/M_PM
68     
69                  DO L = 1, MMAX
70                     LM = FUNLM(L,M)
71                     D_PL = D_P(IJK,L)
72                     M_PL = (PI/6.d0)*(D_PL**3)*RO_S(IJK,L)
73     
74                     MPSUM = M_PM + M_PL
75                     DPSUMo2 = (D_PM+D_PL)/2.d0
76                     NU_PL = ROP_S(IJK,L)/M_PL
77     
78                     ED_common_term = (3.d0/4.d0)*(DPSUMo2*DPSUMo2)*&
79                        (1.d0+C_E)*G_0(IJK,M,L)*NU_PM*NU_PL*(M_PM*&
80                        M_PL/MPSUM)*((M_PM*M_PL)**1.5)
81     
82                     IF (M .eq. L) THEN
83                        Ap_lm = MPSUM/(2.d0)
84                        Dp_lm = M_PL*M_PM/(2.d0*MPSUM)
85                        R1p_lm = 1.d0/( (Ap_lm**1.5)*(Dp_lm**3) )
86                        R3p_lm = 1.d0/( (Ap_lm**1.5)*(Dp_lm**3.5) )
87     
88     ! Dissipation associated with the difference in temperature
89     ! (e.g. interphase transfer term).  For the case case (L=M)
90     ! the term ED_s * (TL-TM) cancels. Therefore, explicity set
91     ! the term to zero for this case.
92                        ED_ss_ip(IJK,LM) = ZERO
93     
94     ! Dissipation associated with temperature
95                        EDT_s_ip(IJK,M,L) = -ED_common_term* (1.d0-C_E)*&
96                           (M_PL/MPSUM)*(DSQRT(PI)/6.d0)*R1p_lm*&
97                           (Theta_m(IJK,M)**1.5)
98     
99     ! Dissipation associated with divergence of
100     ! velocity of solid phase L: do not explicity include
101     ! terms that cancel when EDvel_sL is summed with EDvel_sM
102                        EDvel_sL = ED_common_term*((1.d0-C_E)*(M_PL/MPSUM)*&
103                           (DPSUMo2*PI/48.d0)*(M_PM*M_PL/MPSUM)*R3p_lm)
104                           EDvel_sL_ip(IJK,M,L) = EDvel_sL*(Theta_m(IJK,L))
105     
106     ! Dissipation associated with divergence of
107     ! velocity of solid phase M: do not explicity include
108     ! terms that cancel when EDvel_sL is summed with EDvel_sM
109     ! commented by sof, no need to re-do computation
110                        !EDvel_sM = ED_common_term*((1.d0-C_E)*(M_PL/MPSUM)*&
111                        !   (DPSUMo2*PI/48.d0)*(M_PM*M_PL/MPSUM)*R3p_lm)
112                        !EDvel_sM_ip(IJK,M,L) = EDvel_sM*(Theta_m(IJK,M)) ! M=L
113     
114                        EDvel_sM_ip(IJK,M,L) =  EDvel_sL_ip(IJK,M,L)
115     
116                     ELSE
117     
118                        Ap_lm = (M_PM*Theta_m(IJK,L)+M_PL*Theta_m(IJK,M))/&
119                           2.d0
120                        Bp_lm = (M_PM*M_PL*(Theta_m(IJK,L)-&
121                           Theta_m(IJK,M) ))/(2.d0*MPSUM)
122                        Dp_lm = (M_PL*M_PM*(M_PM*Theta_m(IJK,M)+M_PL*&
123                           Theta_m(IJK,L) ))/(2.d0*MPSUM*MPSUM)
124     
125                        R1p_lm = (1.d0/((Ap_lm**1.5)*(Dp_lm**3)))+ &
126                           ((9.d0*Bp_lm*Bp_lm)/(Ap_lm**2.5 * Dp_lm**4))+&
127                           ((30.d0*Bp_lm**4)/(2.d0*Ap_lm**3.5 * Dp_lm**5))
128     
129                        R3p_lm = (1.d0/((Ap_lm**1.5)*(Dp_lm**3.5)))+&
130                           ((21.d0*Bp_lm*Bp_lm)/(2.d0 * Ap_lm**2.5 * Dp_lm**4.5))+&
131                           ((315.d0*Bp_lm**4)/(8.d0 * Ap_lm**3.5 *Dp_lm**5.5))
132     
133                        R4p_lm = (3.d0/( Ap_lm**2.5 * Dp_lm**3.5))+&
134                           ((35.d0*Bp_lm*Bp_lm)/(2.d0 * Ap_lm**3.5 * Dp_lm**4.5))+&
135                           ((441.d0*Bp_lm**4)/(8.d0 * Ap_lm**4.5 * Dp_lm**5.5))
136     
137                        R5p_lm = (1.d0/(Ap_lm**2.5 * Dp_lm**3))+ &
138                           ((5.d0*Bp_lm*Bp_lm)/(Ap_lm**3.5 * Dp_lm**4))+&
139                           ((14.d0*Bp_lm**4)/( Ap_lm**4.5 * Dp_lm**5))
140     
141                        R10p_lm = (1.d0/(Ap_lm**2.5 * Dp_lm**2.5))+&
142                           ((25.d0*Bp_lm*Bp_lm)/(2.d0* Ap_lm**3.5 * Dp_lm**3.5))+&
143                           ((1225.d0*Bp_lm**4)/(24.d0* Ap_lm**4.5 * Dp_lm**4.5))
144     
145     ! Dissipation associated with the difference in temperature
146     ! (e.g. interphase transfer term). to solved using PEA
147                        ED_ss_ip(IJK,LM) = ED_common_term*DSQRT(PI)*(M_PM*M_PL/&
148                           (2.d0*MPSUM))*R5p_lm*( (Theta_M(IJK,M)*Theta_M(IJK,L))**3 )
149     
150     ! Dissipation associated with temperature
151                        EDT_s_ip(IJK,M,L) = -ED_common_term* (1.d0-C_E)*(M_PL/MPSUM)*&
152                           (DSQRT(PI)/6.d0)*R1p_lm*( (Theta_m(IJK,M)*Theta_m(IJK,L))**3 )
153     
154     ! Dissipation associated with divergence of
155     ! velocity of solid phase L
156                        EDvel_sL = ED_common_term*( ((3.d0*DPSUMo2*PI/40.d0)*M_PL*&
157                           R10p_lm)+( (DPSUMo2*PI/4.d0)*(M_PM*M_PL/MPSUM)*Bp_lm*&
158                           R4p_lm)-( (1.d0-C_E)*(M_PL/MPSUM)*(DPSUMo2*PI/16.d0)*&
159                           M_PL*Bp_lm*R4p_lm)+( (1.d0-C_E)*(M_PL/MPSUM)*&
160                           (DPSUMo2*PI/48.d0)*(M_PM*M_PL/MPSUM)*R3p_lm))
161     
162                        EDvel_sL_ip(IJK,M,L) = EDvel_sL*( Theta_m(IJK,M)**3.5 *&
163                           Theta_m(IJK,L)**2.5 )
164     
165     ! Dissipation associated with divergence of
166     ! velocity of solid phase M
167                        EDvel_sM = ED_common_term*( (-(3.d0*DPSUMo2*PI/40.d0)*M_PM*&
168                           R10p_lm)+( (DPSUMo2*PI/4.d0)*(M_PM*M_PL/MPSUM)*Bp_lm*&
169                           R4p_lm)+( (1.d0-C_E)*(M_PL/MPSUM)*(DPSUMo2*PI/16.d0)*&
170                            M_PM*Bp_lm*R4p_lm)+( (1.d0-C_E)*(M_PL/MPSUM)*&
171                           (DPSUMo2*PI/48.d0)*(M_PM*M_PL/MPSUM)*R3p_lm))
172     
173                        EDvel_sM_ip(IJK,M,L) = EDvel_sM*( Theta_m(IJK,M)**2.5 *&
174                           Theta_m(IJK,L)**3.5 )
175     
176                     ENDIF   ! end if/else (m.eq.l)
177     
178                  ENDDO   ! end do (l=1,mmax)
179     
180              ENDIF   ! end if(fluid_at)
181           ENDDO   ! end do ijk
182     
183           RETURN
184           END SUBROUTINE CALC_IA_ENERGY_DISSIPATION_SS
185     
186     
187     
188     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
189     !                                                                      C
190     !  Subroutine: CALC_GD_99_ENERGY_DISSIPATION_SS                        C
191     !                                                                      C
192     !  Purpose: Implement kinetic theory of Garzo & Dufty (1999)           C
193     !     for calculation of source terms in granular energy equation      C
194     !                                                                      C
195     !  Author: Janine E. Galvin                                            C
196     !                                                                      C
197     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
198     
199           SUBROUTINE CALC_GD_99_ENERGY_DISSIPATION_SS(M, IER)
200     
201     !-----------------------------------------------
202     ! Modules
203     !-----------------------------------------------
204           USE compar
205           USE constant
206           USE fldvar
207           USE fun_avg
208           USE functions
209           USE geometry
210           USE indices
211           USE param
212           USE param1
213           USE physprop
214           USE rdf
215           USE run
216           USE toleranc
217           use kintheory
218           IMPLICIT NONE
219     !-----------------------------------------------
220     ! Dummy arguments
221     !-----------------------------------------------
222     ! Error index
223           INTEGER, INTENT(INOUT) :: IER
224     ! Solids phase index
225           INTEGER, INTENT(IN) :: M
226     !-----------------------------------------------
227     ! Local variables
228     !-----------------------------------------------
229     ! Indices
230           INTEGER :: IJK, I, J, K
231     ! variables for GD model
232           DOUBLE PRECISION :: press_star, c_star, zeta0_star, nu_gamma_star, &
233                               lambda_num, cd_num, zeta1
234           DOUBLE PRECISION :: D_PM, EP_SM
235           DOUBLE PRECISION :: nu0, Chi
236     !-----------------------------------------------
237     
238           DO IJK = ijkstart3, ijkend3
239               I = I_OF(IJK)
240               J = J_OF(IJK)
241               K = K_OF(IJK)
242     
243               IF ( FLUID_AT(IJK) ) THEN
244     
245     ! Note: k_boltz = M_PM
246     
247     ! local aliases
248                  Chi = G_0(IJK,M,M)
249                  EP_SM = EP_s(IJK,M)
250                  D_PM = D_P(IJK,M)
251     
252     !            nu0=p_k/eta0
253                  nu0 = (96.d0/5.d0)*(EP_SM/D_PM)*DSQRT(Theta_m(IJK,M)/PI)
254     
255                  press_star = 1.d0 + 2.d0*(1.d0+C_E)*EP_SM*Chi
256     
257                  c_star = 32.0d0*(1.0d0 - C_E)*(1.d0 - 2.0d0*C_E*C_E) &
258                     / (81.d0 - 17.d0*C_E + 30.d0*C_E*C_E*(1.0d0-C_E))
259     
260                  zeta0_star = (5.d0/12.d0)*Chi*(1.d0 - C_E*C_E) &
261                     * (1.d0 + (3.d0/32.d0)*c_star)
262     
263                  nu_gamma_star = ((1.d0+C_E)/48.d0)*Chi*(128.d0-96.d0*C_E + &
264                     15.d0*C_E*C_E-15.d0*C_E*C_E*C_E+ (c_star/64.d0) * (15.d0* &
265                     C_E*C_E*C_E-15.d0*C_E*C_E+498.d0*C_E-434.d0))
266     
267                  lambda_num = (3.d0/8.d0)*( (1.d0-C_E)*(5.d0*C_E*C_E+4.d0*C_E-1.d0)+ &
268                     (c_star/12.d0)*(159.d0*C_E+3.d0*C_E*C_E-19.d0*C_E- &
269                     15.d0*C_E*C_E*C_E) ) * (1.d0+C_E)
270     
271     ! does not include factor of 1/nu0.  the factor 1/nu0 in cD will cancel
272     ! with the factor of nu0 used to obtain zeta1 from zeta1_star
273                  cd_num = ( (4.d0/15.d0)*lambda_num*EP_SM*Chi + &
274                     (press_star-1.d0)*(2.d0/3.d0-C_E)*c_star ) / &
275                     ( 0.5d0*zeta0_star+nu_gamma_star + (5.d0*c_star/64.d0) * &
276                     (1.d0+(3.d0*c_star/64.d0))*Chi*(1.d0-C_E*C_E))
277     
278     ! does not include factor of 1/nu0 in the first term.  the factor 1/nu0
279     ! will cancel with the factor of nu0 used to obtain zeta1 from zeta1_star
280     ! zeta1 = nu0*zeta1_star
281                  zeta1 = -(1.d0-C_E)*(press_star-1.d0) + (5.d0/32.d0) * &
282                     (1.d0-C_E*C_E)*(1.d0+(3.d0*c_star/64.d0))*Chi*cd_num
283     
284     ! in the energy equation this term is multiplied by 3/2*n*kboltz*T
285     ! leave multiplication of theta for source routine
286                  EDvel_sM_ip(IJK,M,M) = (3.d0/2.d0)*ROP_s(IJK,M)*zeta1
287     
288     ! in the energy equation this term is multiplied by 3/2*n*kboltz*T
289     ! leave multiplication of theta for source routine
290                  EDT_s_ip(IJK,M,M) = (3.d0/2.d0)*ROP_s(IJK,M)*nu0*zeta0_star
291     
292              ENDIF   ! end if (fluid_at)
293           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
294     
295           RETURN
296           END SUBROUTINE CALC_GD_99_ENERGY_DISSIPATION_SS
297     
298     
299     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
300     !                                                                      C
301     !  Subroutine: CALC_GTSH_ENERGY_DISSIPATION_SS                         C
302     !                                                                      C
303     !  Purpose: Implement kinetic theory of Garzo, Tenneti, Subramaniam    C
304     !     Hrenya (2012) for calculation of source terms in granular        C
305     !     energy equation                                                  C
306     !                                                                      C
307     !  Author: Sofiane Benyahia                                            C
308     !                                                                      C
309     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
310     
311           SUBROUTINE CALC_GTSH_ENERGY_DISSIPATION_SS(M)
312     
313     !-----------------------------------------------
314     ! Modules
315     !-----------------------------------------------
316           USE compar
317           USE constant
318           USE fldvar
319           USE fun_avg
320           USE functions
321           USE geometry
322           USE indices
323           USE param
324           USE param1
325           USE physprop
326           USE rdf
327           USE run
328           USE toleranc
329           use kintheory
330           IMPLICIT NONE
331     !-----------------------------------------------
332     ! Dummy arguments
333     !-----------------------------------------------
334     ! Solids phase index
335           INTEGER, INTENT(IN) :: M
336     !-----------------------------------------------
337     ! Local variables
338     !-----------------------------------------------
339     ! Indices
340           INTEGER :: IJK, I, J, K, IMJK, IJMK, IJKM
341     ! variables
342           DOUBLE PRECISION :: D_PM, EP_SM, V_p, N_p, M_p
343           DOUBLE PRECISION :: nu0, Chi
344           DOUBLE PRECISION :: UGC, VGC, WGC, USCM, VSCM, WSCM, VREL
345           DOUBLE PRECISION :: Re_m, Re_T
346           DOUBLE PRECISION :: zeta_star, mu2_0, mu4_0, mu4_1
347           DOUBLE PRECISION :: omega, nu_j, rho_10, rho_11
348     
349     !-----------------------------------------------
350     ! Function subroutines
351     !-----------------------------------------------
352           DOUBLE PRECISION S_star
353           DOUBLE PRECISION G_gtsh
354     !-----------------------------------------------
355     !
356           DO IJK = ijkstart3, ijkend3
357               I = I_OF(IJK)
358               J = J_OF(IJK)
359               K = K_OF(IJK)
360               IMJK = IM_OF(IJK)
361               IJMK = JM_OF(IJK)
362               IJKM = KM_OF(IJK)
363     
364               IF ( FLUID_AT(IJK) ) THEN
365     
366     ! Local aliases
367                  Chi = G_0(IJK,M,M)
368                  EP_SM = EP_s(IJK,M)
369                  D_PM = D_P(IJK,M)
370                  V_p = pi*D_PM**3/6d0
371                  n_p = EP_SM/V_p
372                  M_p = V_p * ro_s(ijk,m)
373     
374                  nu0 = (96.d0/5.d0)*(EP_SM/D_PM)*DSQRT(Theta_m(IJK,M)/PI)
375     
376     
377     ! First calculate the Re number (Re_m, Re_T) in eq. 3.1 in GTSH theory
378                 UGC = AVG_X_E(U_G(IMJK),U_G(IJK),I)
379                 VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
380                 WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
381                 USCM = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
382                 VSCM = AVG_Y_N(V_S(IJMK,M),V_S(IJK,M))
383                 WSCM = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M))
384                 VREL = DSQRT((UGC - USCM)**2 + (VGC - VSCM)**2 + (WGC - WSCM)**2)
385     
386                 Re_m = D_PM*VREL*ROP_g(ijk)/Mu_g(ijk)  ! rop_g = ro_g * (1-phi)
387                 Re_T = ro_g(ijk)*D_PM*dsqrt(theta_m(ijk,m)) / mu_g(ijk)
388     
389     ! Now calculate and store xsi (used in many spots) in eq. 8.2 in
390     ! GTSH theory.
391                 xsi_gtsh(ijk) = one/6d0*D_PM*VREL**2*(3d0*pi*mu_g(ijk)*&
392                    D_PM/M_p)**2 / dsqrt(pi*theta_m(ijk,m)) * &
393                    S_star(EP_SM, Chi)
394     
395     ! eq. (6.22) GTSH theory
396                 mu2_0 = dsqrt(2d0*pi) * Chi * (one-C_E**2)
397     
398     ! eq. (6.23) GTSH theory
399                 mu4_0 = (4.5d0+C_E**2) * mu2_0
400     
401     ! eq. (6.24) GTSH theory
402     !            mu4_1 = (6.46875d0+0.3125d0*C_E**2 + 2d0/(one-C_E)) * mu2_0
403     ! this is done to avoid /0 in case c_e = 1.0
404                 mu4_1 = (6.46875d0+0.9375d0*C_E**2)*mu2_0 + &
405                    2d0*dsqrt(2d0*pi)* Chi*(one+C_E)
406     
407                 A2_gtsh(ijk) = zero ! for EP_SM = zero
408                 if(EP_SM> small_number) then ! avoid singularity
409     ! Now calculate zeta_star in eq. 8.10 in GTSH theory
410                    zeta_star = 4.5d0*dsqrt(2d0*Pi)*&
411                       (ro_g(ijk)/ro_s(ijk,m))**2*Re_m**2 * &
412                       S_star(EP_SM,Chi) / (EP_SM*(one-EP_SM)**2 * Re_T**4)
413     
414     ! Now calculate important parameter A2. This is used in many transport
415     ! coefficients so that it's best to store it in an array instead of
416     ! having many function calls.
417                    A2_gtsh(ijk) = (5d0*mu2_0 - mu4_0) / (mu4_1 - 5d0* &
418                       (19d0/16d0*mu2_0 - 1.5d0*zeta_star))
419                 endif
420     
421     ! Now calculate first term rho_0 in cooling rate rho, eq (6.26) in GTSH
422     ! theory.
423     ! note that theta_(ijk,m) does't contain mass m, so that T/m = theta_m.
424     ! note that edt_s_ip will need to be multiplied by 3/2 rop_s(ijk,m) in
425     ! source_granular_energy to have the same meaning as in GD_99 theory.
426                 EDT_s_ip(ijk,M,M) = 4d0/3d0*dsqrt(pi)*(one-C_E**2)*Chi* &
427                    (one+0.1875d0*A2_gtsh(ijk))*n_p*D_PM**2*dsqrt(theta_m(ijk,m))
428     
429     ! Calculate the second term in cooling rate, eq (7.25) in GTSH theory.
430     ! First calculate eq (7.28, 7.29) of GTSH theory.
431                 omega = (one+C_E)*nu0*((one-C_E**2)*(5d0*C_E-one) - &
432                    A2_gtsh(ijk)/ 6d0 * (15d0*C_E**3-3d0*C_E**2+81d0*C_E-61d0))
433                 nu_j = (one+C_E)/192d0*Chi*nu0* &
434                    (241d0-177d0*C_E+30d0*C_E**2-30d0*C_E**3)
435     
436     ! Now calculate eq (7.26, 7.27) of GTSH theory. ! corrected by W. Fullmer
437                 rho_10 = 2d0*Chi*EP_SM*(C_E**2-one)
438                 rho_11 = 25d0/1024d0*EP_SM*Chi**2*(one-C_E**2)* &
439                    (one+3d0/128d0*A2_gtsh(ijk)) * (omega/10d0 - &
440                    (one+C_E)*nu0*(one/3d0-C_E)*A2_gtsh(ijk)/2d0) / &
441                    (nu_j + G_gtsh(EP_SM, Chi, IJK, M)/m_p + 1.5d0* &
442                    xsi_gtsh(ijk)/theta_m(ijk,m) -1.5d0*EDT_s_ip(ijk,M,M))
443     
444     ! note that EDvel_sM_ip will later need to be multiplied by 3/2 rop_s(ijk,m).
445                 EDvel_sM_ip(IJK,M,M) = rho_10 + rho_11
446     
447              ENDIF   ! end if (fluid_at)
448           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
449     
450           RETURN
451           END SUBROUTINE CALC_GTSH_ENERGY_DISSIPATION_SS
452     
453     
454     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
455     !                                                                      C
456     !  Functions: G_tsh, K_phi, R_d, S_star                                C
457     !                                                                      C
458     !  Purpose: Implement kinetic theory of Garzo, Tenneti, Subramaniam    C
459     !     Hrenya (2012) for calculation of source terms in granular        C
460     !     energy equation                                                  C
461     !                                                                      C
462     !  Author: Sofiane Benyahia                                            C
463     !                                                                      C
464     !  Comments:                                                           C
465     !     these functions are needed only for gtsh theory                  C
466     !     Function gamma, eq. (8.1) in GTSH theory                         C
467     !                                                                      C
468     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
469     
470           DOUBLE PRECISION FUNCTION G_gtsh (EP_SM, Chi, IJK, M)
471     
472     !-----------------------------------------------
473     ! Modules
474     !-----------------------------------------------
475           USE param
476           USE param1
477           USE physprop
478           USE fldvar
479           USE constant
480           IMPLICIT NONE
481     !-----------------------------------------------
482     ! Dummy arguments
483     !-----------------------------------------------
484     ! solids volume fraction of phase m at ijk
485           DOUBLE PRECISION, INTENT(IN) :: EP_SM
486     ! radial distribution function of phase M at ijk
487           DOUBLE PRECISION, INTENT(IN) :: Chi
488     ! Index
489           INTEGER, INTENT(IN) :: IJK
490     ! Solids phase index. Note M should be equal to 1 since theory
491     ! valid for only mmax = 1.
492           INTEGER, INTENT(IN) :: M
493     !-----------------------------------------------
494     ! Local variables
495     !-----------------------------------------------
496           DOUBLE PRECISION :: Re_T
497           DOUBLE PRECISION :: Rdiss, RdissP
498     !-----------------------------------------------
499     ! Functions
500     !-----------------------------------------------
501           DOUBLE PRECISION :: K_phi
502     !-----------------------------------------------
503     
504           if(EP_SM <= 0.1d0) then
505              RdissP = one+3d0*dsqrt(EP_SM/2d0)
506           else
507              RdissP = &
508              one + 3d0*dsqrt(EP_SM/2d0) + 135d0/64d0*EP_SM*dlog(EP_SM) + &
509              11.26d0*EP_SM*(one-5.1d0*EP_SM+16.57d0*EP_SM**2-21.77d0*    &
510              EP_SM**3) - EP_SM*Chi*dlog(epM)
511           endif
512     
513           Re_T = ro_g(ijk)*d_p(ijk,m)*dsqrt(theta_m(ijk,m)) / mu_g(ijk)
514     
515           Rdiss = RdissP + Re_T * K_phi(EP_SM)
516     
517           G_gtsh = 3d0*Pi*mu_g(ijk)*d_p(ijk,m)*Rdiss  ! eq. (8.1) in GTSH theory
518     
519           RETURN
520           END FUNCTION G_gtsh
521     
522     
523     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
524     !                                                                      C
525     !  Comments:                                                           C
526     !     Function K(phi), eq. (2.7) in Yi/Zenk/Mitrano/Hrenya JFM (2013)  C
527     !                                                                      C
528     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
529           DOUBLE PRECISION FUNCTION K_phi (phi)
530           IMPLICIT NONE
531     !-----------------------------------------------
532     ! Dummy Arguments
533     !-----------------------------------------------
534     ! solids volume fraction
535           DOUBLE PRECISION, INTENT(IN) :: phi
536     !-----------------------------------------------
537     
538           K_phi = (0.096d0 + 0.142d0*phi**0.212d0) / (1d0-phi)**4.454d0
539           K_phi = 0.0d0 ! set to zero for compatibility with GTSH JFM (2012)
540           RETURN
541           END FUNCTION K_phi
542     
543     
544     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
545     !                                                                      C
546     !  Comments:                                                           C
547     !     Function R_d, eq. (8.6) in GTSH theory                           C
548     !                                                                      C
549     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
550           DOUBLE PRECISION FUNCTION R_d (phi)
551           IMPLICIT NONE
552     !-----------------------------------------------
553     ! Dummy arguments
554     !-----------------------------------------------
555     ! solids volume fraction
556           DOUBLE PRECISION, INTENT(IN) :: phi
557     !-----------------------------------------------
558     
559           R_d = 1.0d0  ! this avoids singularity at phi = 0.0
560           if((phi > 1d-15) .and. (phi <= 0.4d0)) then
561             R_d = (1d0+3d0*dsqrt(phi/2d0)+135d0/64d0*phi*dlog(phi)+17.14d0*phi) / &
562                   (1d0+0.681d0*phi-8.48*phi**2+8.16d0*phi**3)
563           elseif(phi > 0.4d0) then
564             R_d = 10d0*phi/(1d0-phi)**3 + 0.7d0
565           endif
566           RETURN
567           END FUNCTION R_d
568     
569     
570     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
571     !                                                                      C
572     !  Comments:                                                           C
573     !     Function S_Star(phi), eq. (8.3, 8.5) in GTSH theory              C
574     !                                                                      C
575     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
576           DOUBLE PRECISION FUNCTION S_star (phi, Chi)
577           IMPLICIT NONE
578     !-----------------------------------------------
579     ! Dummy arguments
580     !-----------------------------------------------
581     ! solids volume fraction
582           DOUBLE PRECISION, INTENT(IN) :: phi
583     ! radial distribution function
584           DOUBLE PRECISION, INTENT(IN) :: Chi
585     !-----------------------------------------------
586     ! Functions
587     !-----------------------------------------------
588           DOUBLE PRECISION R_d
589     !-----------------------------------------------
590     
591           S_star = 1.0d0
592           if(phi >= 0.1d0) &
593             S_star = R_d(phi)**2/(Chi*(1d0+3.5d0*dsqrt(phi)+5.9*phi))
594           RETURN
595           END FUNCTION S_Star
596