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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOURCE_GRANULAR_ENERGY                                  C
4     !  Purpose: Calculate the source terms in the granular energy equation C
5     !  Note: The center coefficient and source vector are negative.        C
6     !  The off-diagonal coefficients are positive.                         C
7     !                                                                      C
8     !  Author: Kapil Agrawal, Princeton University        Date: 04-FEB-98  C
9     !  Reviewer: M. Syamlal                               Date:            C
10     !                                                                      C
11     !  Revision Number:1                                                   C
12     !  Purpose: Add Simonin and Ahmadi models                              C
13     !  Author: Sofiane Benyahia, Fluent Inc.               Date: 02-01-05  C
14     !                                                                      C
15     !  Literature/Document References:                                     C
16     !     Lun, C.K.K., S.B. Savage, D.J. Jeffrey, and N. Chepurniy,        C
17     !        Kinetic theories for granular flow - inelastic particles in   C
18     !        Couette-flow and slightly inelastic particles in a general    C
19     !        flow field. Journal of Fluid Mechanics, 1984. 140(MAR):       C
20     !        p. 223-256                                                    C
21     !     Gidaspow, D., Multiphase flow and fluidziation, 1994, Academic   C
22     !        Press Inc., California, Chapter 9.                            C
23     !     Koch, D. L., and Sangani, A. S., Particle pressure and marginal  C
24     !        stability limits for a homogeneous monodisperse gas-fluidized C
25     !        bed: kinetic theory and numerical simulations, Journal of     C
26     !        Fluid Mechanics, 1999, 400, 229-263.                          C
27     !                                                                      C
28     !     Simonin, O., 1996. Combustion and turbulence in two-phase flows, C
29     !        Von Karman institute for fluid dynamics, lecture series,      C
30     !        1996-02                                                       C
31     !     Balzer, G., Simonin, O., Boelle, A., and Lavieville, J., 1996,   C
32     !        A unifying modelling approach for the numerical prediction    C
33     !        of dilute and dense gas-solid two phase flow. CFB5, 5th int.  C
34     !        conf. on circulating fluidized beds, Beijing, China.          C
35     !     Cao, J. and Ahmadi, G., 1995, Gas-particle two-phase turbulent   C
36     !        flow in a vertical duct. Int. J. Multiphase Flow, vol. 21,    C
37     !        No. 6, pp. 1203-1228.                                         C
38     !                                                                      C
39     !                                                                      C
40     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
41     
42           SUBROUTINE SOURCE_GRANULAR_ENERGY(SOURCELHS, &
43                         SOURCERHS, IJK, M)
44     
45     !-----------------------------------------------
46     ! Modules
47     !-----------------------------------------------
48           USE compar
49           USE constant
50           USE drag
51           USE fldvar
52           USE fun_avg
53           USE functions
54           USE geometry
55           USE indices
56           USE kintheory
57           USE mms
58           USE parallel
59           USE param
60           USE param1
61           USE physprop
62           USE rdf
63           USE run
64           USE solids_pressure
65           USE toleranc
66           USE trace
67           USE turb
68           USE visc_g
69           USE visc_s
70     
71           IMPLICIT NONE
72     !-----------------------------------------------
73     ! Dummy Arguments
74     !-----------------------------------------------
75     ! Source terms to be kept on rhs (source vector)
76           DOUBLE PRECISION, INTENT(INOUT) :: sourcerhs
77     ! Source terms to be kept on lhs (center coefficient)
78           DOUBLE PRECISION, INTENT(INOUT) :: sourcelhs
79     ! Solid phase index
80           INTEGER, INTENT(IN) :: M
81     ! Indices
82           INTEGER, INTENT(IN) :: IJK
83     !-----------------------------------------------
84     ! Local variables
85     !-----------------------------------------------
86     ! Indices
87           INTEGER :: I, J, K, IM, JM, KM, IMJK, IJMK, IJKM, &
88                      IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
89     ! Solid phase index
90           INTEGER :: MM
91     ! Particle relaxation time
92           DOUBLE PRECISION :: Tau_12_st
93     ! Sum of eps*G_0
94           DOUBLE PRECISION :: SUM_EpsGo
95     ! Slip velocity
96           DOUBLE PRECISION :: VSLIP
97     ! coefficients in heat flux term
98           DOUBLE PRECISION :: Knu_e, Knu_w, Knu_n, Knu_s, &
99                               Knu_t, Knu_b
100     ! Source terms to be kept on rhs
101           DOUBLE PRECISION :: S1_rhs, S2_rhs, S3_rhs, &
102                               S5a_rhs, S5b_rhs, S5c_rhs, S5_rhs, &
103                               S7_rhs
104     ! Source terms to be kept on lhs
105           DOUBLE PRECISION :: S1_lhs, S3_lhs, S4_lhs, S6_lhs
106     !-----------------------------------------------
107     
108           I = I_OF(IJK)
109           J = J_OF(IJK)
110           K = K_OF(IJK)
111           IM = Im1(I)
112           JM = Jm1(J)
113           KM = Km1(K)
114           IMJK = IM_OF(IJK)
115           IJMK = JM_OF(IJK)
116           IJKM = KM_OF(IJK)
117           IJKE = EAST_OF(IJK)
118           IJKW = WEST_OF(IJK)
119           IJKN = NORTH_OF(IJK)
120           IJKS = SOUTH_OF(IJK)
121           IJKT = TOP_OF(IJK)
122           IJKB = BOTTOM_OF(IJK)
123     
124     ! initialize summation variables
125           S1_rhs = ZERO
126           S1_lhs = ZERO
127           S2_rhs = ZERO
128           S3_rhs = ZERO
129           S3_lhs = ZERO
130           S4_lhs = ZERO
131           S6_lhs = ZERO
132           S7_rhs = ZERO
133     
134     ! Changes needed for multitype particles, sof June 16 2005
135     ! Sum of eps*G_0 is used instead of Eps*G_0
136           SUM_EpsGo = ZERO
137           DO MM = 1, MMAX
138              SUM_EpsGo =  SUM_EpsGo+EP_s(IJK,MM)*G_0(IJK,MM,MM)
139           ENDDO
140     
141     ! Production by shear: (S:grad(vi))
142     ! Pi_s*tr(Di)     (Lun et al. 1984)
143           S1_rhs = P_S_C(IJK,M)*ZMAX(( -TRD_S_C(IJK,M) ))
144           S1_lhs = P_S_C(IJK,M)*ZMAX(( TRD_S_C(IJK,M) ))
145     
146     ! Production by shear: (S:grad(vi))
147     ! Mu_s*tr(Di^2)     (Lun et al. 1984)
148           S2_rhs = 2.d0*MU_S_C(IJK,M)*TRD_S2(IJK,M)
149     
150     ! Production by shear: (S:grad(vi))
151     ! Lambda_s*tr(Di)^2     (Lun et al. 1984)
152           S3_rhs = (TRD_S_C(IJK,M)**2)*ZMAX( LAMBDA_S_C(IJK,M) )
153           S3_lhs = (TRD_S_C(IJK,M)**2)*ZMAX( -LAMBDA_S_C(IJK,M) )
154     
155     ! Energy dissipation by collisions
156           S4_lhs = (48.d0/DSQRT(PI))*ETA*(ONE-ETA)*ROP_S(IJK,M)*&
157              SUM_EpsGo*DSQRT(THETA_M(IJK,M))/D_P(IJK,M)
158     
159     ! Need to revisit this part about granular energy MMS source terms.
160           IF (USE_MMS) S4_lhs = ZERO
161     
162     ! Energy dissipation by viscous dampening
163     ! Gidaspow (1994)  : addition due to role of interstitial fluid
164           IF(KT_TYPE_ENUM == SIMONIN_1996 .OR. &
165              KT_TYPE_ENUM == AHMADI_1995) THEN
166              S6_lhs = 3.d0 * F_GS(IJK,M)
167           ELSE IF(SWITCH > ZERO .AND. RO_g0 /= ZERO) THEN
168              S6_lhs = SWITCH * 3.d0 * F_GS(IJK,M)
169           ENDIF
170     
171     ! Energy production due to gas-particle slip
172     ! Koch & Sangani (1999) : addition due to role of interstitial fluid
173           IF(KT_TYPE_ENUM == SIMONIN_1996) THEN
174              S7_rhs = F_GS(IJK,M)*K_12(IJK)
175           ELSEIF(KT_TYPE_ENUM == AHMADI_1995) THEN
176     ! note specific reference to F_GS of solids phase 1!
177              IF(Ep_s(IJK,M) > DIL_EP_S .AND. F_GS(IJK,1) > small_number) THEN
178                 Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,1)
179                 S7_rhs = 2.d0*F_GS(IJK,M) * (ONE/(ONE+Tau_12_st/  &
180                    (Tau_1(ijk)+small_number)))*K_Turb_G(IJK)
181              ELSE
182                 S7_rhs = ZERO
183              ENDIF
184           ELSEIF(SWITCH > ZERO .AND. RO_g0 /= ZERO) THEN
185              VSLIP = DSQRT( (U_S(IJK,M)-U_G(IJK))**2 + &
186                 (V_S(IJK,M)-V_G(IJK))**2 + (W_S(IJK,M)-W_G(IJK))**2 )
187              S7_rhs = SWITCH*81.d0*EP_S(IJK,M)*(MU_G(IJK)*VSLIP)**2 /&
188                 (G_0(IJK,M,M)*D_P(IJK,M)**3 * RO_S(IJK,M)*&
189                 DSQRT(PI)*DSQRT( THETA_M(IJK,M)+SMALL_NUMBER ) )
190     ! Need to revisit this part about granular energy MMS source terms.
191              IF (USE_MMS) S7_rhs = ZERO
192     
193           ENDIF
194     
195     
196     ! The following lines are essentially commented out since Kphi_s has
197     ! been set to zero in subroutine calc_mu_s.  To activate the feature
198     ! activate the following lines and the lines in calc_mu_s.
199     ! Part of Heat Flux: div (q)
200     ! Kphi_s*grad(eps)     (Lun et al. 1984)
201           IF (.FALSE.) THEN
202              Knu_e = AVG_X_H(Kphi_s(IJK,M), Kphi_s(IJKE,M),I)
203              Knu_w = AVG_X_H(Kphi_s(IJKW,M),Kphi_s(IJK,M),IM)
204              Knu_n = AVG_Y_H(Kphi_s(IJK,M), Kphi_s(IJKN,M),J)
205              Knu_s = AVG_Y_H(Kphi_s(IJKS,M),Kphi_s(IJK,M),JM)
206              Knu_t = AVG_Z_H(Kphi_s(IJK,M), Kphi_s(IJKT,M),K)
207              Knu_b = AVG_Z_H(Kphi_s(IJKB,M),Kphi_s(IJK,M),KM)
208     
209              S5a_rhs = Knu_e*(EP_s(IJKE,M)-EP_s(IJK,M))*&
210                 oDX_E(I)*AYZ(IJK) - &
211                 Knu_w*(EP_s(IJK,M)-EP_s(IJKW,M))*&
212                 oDX_E(IM)*AYZ(IMJK)
213              S5b_rhs = Knu_n*(EP_s(IJKN,M)-EP_s(IJK,M))*&
214                 oDY_N(J)*AXZ(IJK) - &
215                 Knu_s*(EP_s(IJK,M)-EP_s(IJKS,M))*&
216                 oDY_N(JM)*AXZ(IJMK)
217              S5c_rhs = Knu_t*(EP_s(IJKT,M)-EP_s(IJK,M))*&
218                 oX(I)*oDZ_T(K)*AXY(IJK) - &
219                 Knu_b*(EP_s(IJK,M)-EP_s(IJKB,M))*&
220                 oX(I)*oDZ_T(KM)*AXY(IJKM)
221              S5_rhs = S5a_rhs + S5b_rhs + S5c_rhs
222           ENDIF
223           S5_rhs = ZERO
224     
225     
226           SOURCERHS = (S1_rhs + S2_rhs + S3_rhs + S7_rhs)*VOL(IJK) + &
227                S5_rhs
228     
229           SOURCELHS = ( ((S1_lhs + S3_lhs)/(THETA_M(IJK,M)+SMALL_NUMBER)) + &
230               S4_lhs + S6_lhs) * VOL(IJK)
231     
232           RETURN
233           END SUBROUTINE SOURCE_GRANULAR_ENERGY
234     
235     
236     
237     
238     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
239     !                                                                      C
240     !  Subroutine: SOURCE_IA_GRANULAR_ENERGY                               C
241     !  Purpose: Calculate the source terms of the Mth solids phase         C
242     !           granular energy equation                                   C
243     !                                                                      C
244     !  Author: Janine E. Galvin, Univeristy of Colorado                    C
245     !                                                                      C
246     !  Literature/Document References:                                     C
247     !    Iddir, Y.H., Modeling of the multiphase mixture of particles      C
248     !       using the kinetic theory approach, PhD Thesis, Illinois        C
249     !       Institute of Technology, Chicago, Illinois, 2004               C
250     !    Iddir, Y.H., & H. Arastoopour, Modeling of multitype particle     C
251     !       flow using the kinetic theory approach, AIChE J., Vol 51,      C
252     !       No 6, June 2005                                                C
253     !                                                                      C
254     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
255     
256           SUBROUTINE SOURCE_GRANULAR_ENERGY_IA(SOURCELHS, &
257                         SOURCERHS, IJK, M)
258     
259     !-----------------------------------------------
260     ! Modules
261     !-----------------------------------------------
262           USE param
263           USE param1
264           USE parallel
265           USE physprop
266           USE run
267           USE drag
268           USE geometry
269           USE fldvar
270           USE visc_g
271           USE visc_s
272           USE trace
273           USE turb
274           USE indices
275           USE constant
276           USE toleranc
277           USE residual
278           use kintheory
279           USE compar
280           USE fun_avg
281           USE functions
282           IMPLICIT NONE
283     !-----------------------------------------------
284     ! Dummy Arguments
285     !-----------------------------------------------
286     ! Source terms to be kept on rhs (source vector)
287           DOUBLE PRECISION, INTENT(INOUT) :: sourcerhs
288     ! Source terms to be kept on lhs (center coefficient)
289           DOUBLE PRECISION, INTENT(INOUT) :: sourcelhs
290     ! Solid phase index
291           INTEGER, INTENT(IN) :: M
292     ! Indices
293           INTEGER, INTENT(IN) :: IJK
294     !-----------------------------------------------
295     ! Local variables
296     !-----------------------------------------------
297     ! Indices
298           INTEGER :: I, J, K, IM, JM, KM, IMJK, IJMK, IJKM,&
299                      IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
300     ! phase index
301           INTEGER :: L, LM
302     ! velocities
303           DOUBLE PRECISION :: UsM_e, UsM_w, VsM_n, VsM_s, WsM_t, WsM_b,&
304                               UsL_e, UsL_w, VsL_n, VsL_s, WsL_t, WsL_b,&
305                               UsM_p, VsM_p, WsM_P, UsL_p, VsL_p, WsL_p
306     ! number densities
307           DOUBLE PRECISION :: NU_PM_E, NU_PM_W, NU_PM_N, NU_PM_S, NU_PM_T,&
308                               NU_PM_B, NU_PM_p,&
309                               NU_PL_E, NU_PL_W, NU_PL_N, NU_PL_S, NU_PL_T,&
310                               NU_PL_B, NU_PL_p
311     ! temperature of species L
312           DOUBLE PRECISION :: T_PL_E, T_PL_W, T_PL_N, T_PL_S, T_PL_T,&
313                               T_PL_B, T_PL_p
314     ! particle characteristics
315           DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL
316     ! coefficients in heat flux term
317           DOUBLE PRECISION :: Knu_sL_e, Knu_sL_w, Knu_sL_n, Knu_sL_s, Knu_sL_t,&
318                               Knu_sL_b, Knu_sM_e, Knu_sM_w, Knu_sM_n, Knu_sM_s,&
319                               Knu_sM_t, Knu_sM_b,&
320                               Kvel_s_e, Kvel_s_w, Kvel_s_n, Kvel_s_s, Kvel_s_t,&
321                               Kvel_s_b,&
322                               Kth_sL_e, Kth_sL_w, Kth_sL_n, Kth_sL_s, Kth_sL_t,&
323                               Kth_sL_b
324     ! Source terms to be kept on rhs
325           DOUBLE PRECISION :: S10_rhs, S15_rhs, S16_rhs,&
326                               S11_sum_rhs, S12_sum_rhs,&
327                               S13_sum_rhs, S17_sum_rhs, S18_sum_rhs
328     ! Source terms
329           DOUBLE PRECISION :: S14a_sum, S14b_sum, S14c_sum, &
330                               S9_sum, s21a_sum, s21b_sum, s21c_sum
331     ! Source terms to be kept on lhs
332           DOUBLE PRECISION :: S10_lhs, S16_lhs,&
333                               S11_sum_lhs, S12_sum_lhs, S13_sum_lhs,&
334                               S17_sum_lhs, S18_sum_lhs, S20_sum_lhs
335     !-----------------------------------------------
336     
337           I = I_OF(IJK)
338           J = J_OF(IJK)
339           K = K_OF(IJK)
340           IM = Im1(I)
341           JM = Jm1(J)
342           KM = Km1(K)
343           IMJK = IM_OF(IJK)
344           IJMK = JM_OF(IJK)
345           IJKM = KM_OF(IJK)
346           IJKE = EAST_OF(IJK)
347           IJKW = WEST_OF(IJK)
348           IJKN = NORTH_OF(IJK)
349           IJKS = SOUTH_OF(IJK)
350           IJKT = TOP_OF(IJK)
351           IJKB = BOTTOM_OF(IJK)
352     
353     ! initialize summation variables
354           S9_sum = ZERO
355           S14a_sum = ZERO
356           S14b_sum = ZERO
357           S14c_sum = ZERO
358           S21a_sum = ZERO
359           S21b_sum = ZERO
360           S21c_sum = ZERO
361     
362           S11_sum_rhs = ZERO
363           S12_sum_rhs = ZERO
364           S13_sum_rhs = ZERO
365           S17_sum_rhs = ZERO
366           S18_sum_rhs = ZERO
367     
368           S11_sum_lhs = ZERO
369           S12_sum_lhs = ZERO
370           S13_sum_lhs = ZERO
371           S17_sum_lhs = ZERO
372           S18_sum_lhs = ZERO
373           S20_sum_lhs = ZERO
374     
375           UsM_e = U_S(IJK,M)
376           UsM_w = U_S(IMJK,M)
377           VsM_n = V_S(IJK,M)
378           VsM_s = V_S(IJMK,M)
379           WsM_t = W_S(IJK,M)
380           WsM_b = W_S(IJKM,M)
381           UsM_p = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
382           VsM_p = AVG_Y_N(V_S(IJMK,M),V_S(IJK,M) )
383           WsM_p = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M) )
384     
385           D_PM = D_P(IJK,M)
386           M_PM = (Pi/6.d0)*D_PM**3 * RO_S(IJK,M)
387           NU_PM_p = ROP_S(IJK,M)/M_PM
388           NU_PM_E = ROP_S(IJKE,M)/M_PM
389           NU_PM_W = ROP_S(IJKW,M)/M_PM
390           NU_PM_N = ROP_S(IJKN,M)/M_PM
391           NU_PM_S = ROP_S(IJKS,M)/M_PM
392           NU_PM_T = ROP_S(IJKT,M)/M_PM
393           NU_PM_B = ROP_S(IJKB,M)/M_PM
394     
395     ! Production by shear: (S:grad(vi))
396     ! Pi_s*tr(Di)
397           S10_lhs = P_S_C(IJK,M) * ZMAX(TRD_S_C(IJK,M))
398           S10_rhs = P_S_C(IJK,M) * ZMAX(-TRD_S_C(IJK,M))
399     
400     
401     ! Production by shear: (S:grad(vi))
402     ! Mu_s*tr(Di^2)
403           S15_rhs = 2.d0*Mu_s_c(IJK,M)*TRD_S2(IJK,M)
404     
405     ! Production by shear: (S:grad(vi))
406     ! Lambda_s*tr(Di)^2
407           S16_lhs = (TRD_S_C(IJK,M)**2)*ZMAX( -LAMBDA_s_c(IJK,M) )
408           S16_rhs = (TRD_S_C(IJK,M)**2)*ZMAX(  LAMBDA_s_C(IJK,M) )
409     
410           DO L = 1, MMAX
411               D_PL = D_P(IJK,L)
412               M_PL = (Pi/6.d0)*D_PL**3 * RO_S(IJK,L)
413               NU_PL_p = ROP_S(IJK,L)/M_PL
414               NU_PL_E = ROP_S(IJKE,L)/M_PL
415               NU_PL_W = ROP_S(IJKW,L)/M_PL
416               NU_PL_N = ROP_S(IJKN,L)/M_PL
417               NU_PL_S = ROP_S(IJKS,L)/M_PL
418               NU_PL_T = ROP_S(IJKT,L)/M_PL
419               NU_PL_B = ROP_S(IJKB,L)/M_PL
420     
421               T_PL_p = Theta_m(IJK,L)
422               T_PL_E = Theta_m(IJKE,L)
423               T_PL_W = Theta_m(IJKW,L)
424               T_PL_N = Theta_m(IJKN,L)
425               T_PL_S = Theta_m(IJKS,L)
426               T_PL_T = Theta_m(IJKT,L)
427               T_PL_B = Theta_m(IJKB,L)
428     
429               UsL_e = U_S(IJK,L)
430               UsL_w = U_S(IMJK,L)
431               VsL_n = V_S(IJK,L)
432               VsL_s = V_S(IJMK,L)
433               WsL_t = W_S(IJK,L)
434               WsL_b = W_S(IJKM,L)
435               UsL_p = AVG_X_E(U_S(IMJK,L),U_S(IJK,L),I)
436               VsL_p = AVG_Y_N(V_S(IJMK,L),V_S(IJK,L) )
437               WsL_p = AVG_Z_T(W_S(IJKM,L),W_S(IJK,L) )
438     
439     ! Energy dissipation by collisions: Sum(Nip)
440     ! SUM( EDT_s_ip )
441               S20_sum_lhs = S20_sum_lhs + EDT_s_ip(IJK,M,L)
442     
443     ! Energy dissipation by collisions: SUM(Nip)
444     ! SUM( EDvel_sL_ip* div(vp) ) !Modified by sof to include trace of V_s_L
445               S11_sum_lhs = S11_sum_lhs + ZMAX(-EDvel_sL_ip(IJK,M,L)* &
446                    TRD_S_C(IJK,L) ) * VOL(IJK)
447               S11_sum_rhs = S11_sum_rhs + ZMAX( EDvel_sL_ip(IJK,M,L)* &
448                    TRD_S_C(IJK,L) ) * VOL(IJK)
449     
450     ! Energy dissipation by collisions: Sum(Nip)
451     ! SUM( EDvel_sM_ip* div(vi) ) !Modified by sof to include trace of V_s_M
452               S12_sum_lhs = S12_sum_lhs + ZMAX(-EDvel_sM_ip(IJK,M,L)*&
453                    TRD_S_C(IJK,M) ) * VOL(IJK)
454               S12_sum_rhs = S12_sum_rhs + ZMAX( EDvel_sM_ip(IJK,M,L)*&
455                    TRD_S_C(IJK,M) ) * VOL(IJK)
456     
457               IF (M .NE. L) THEN
458                    LM = FUNLM(L,M)
459     
460     ! Production by shear: (S:grad(vi))
461     ! SUM(2*Mu_sL_ip*tr(Dk*Di) )
462                    S17_sum_lhs = S17_sum_lhs + 2.d0*MU_sL_ip(IJK,M,L)*&
463                         ZMAX( - TRD_s2_ip(IJK,M,L) )
464                    S17_sum_rhs = S17_sum_rhs + 2.d0*MU_sL_ip(IJK,M,L)*&
465                         ZMAX( TRD_s2_ip(IJK,M,L) )
466     
467     ! These two terms can be treated explicitly here by uncommenting the following
468     ! two lines. They are currently treated by PEA algorithm in solve_granular_energy
469     !
470     !               S13_sum_lhs = S13_sum_lhs + ED_ss_ip(IJK,LM)*Theta_m(IJK,M)
471     !               S13_sum_rhs = S13_sum_rhs + ED_ss_ip(IJK,LM)*Theta_m(IJK,L)
472     !
473     !
474     ! Production by shear: (S:grad(vi))
475     ! SUM( (Xi_sL_ip-(2/3)*Mu_sL_ip)*tr(Dk)tr(Di) )
476                    S18_sum_lhs = S18_sum_lhs + ZMAX(-1.d0* &
477                         (Xi_sL_ip(IJK,M,L)-(2.d0/3.d0)*Mu_sL_ip(IJK,M,L))*&
478                         TRD_S_C(IJK,M)*TRD_S_C(IJK,L) )
479                    S18_sum_rhs = S18_sum_rhs + ZMAX(&
480                         (Xi_sL_ip(IJK,M,L)-(2.d0/3.d0)*Mu_sL_ip(IJK,M,L))*&
481                         TRD_S_C(IJK,M)*TRD_S_C(IJK,L) )
482     
483     ! Part of Heat Flux: div (q)
484     ! Kth_sL_ip*[grad(Tp)]
485     ! Note for L=M S21 terms cancel with similar term arising from grad(Ti)
486                    Kth_sL_e = AVG_X_S(Kth_sL_ip(IJK,M,L), Kth_sL_ip(IJKE,M,L),I)
487                    Kth_sL_w = AVG_X_S(Kth_sL_ip(IJKW,M,L),Kth_sL_ip(IJK,M,L), IM)
488                    Kth_sL_n = AVG_Y_S(Kth_sL_ip(IJK,M,L), Kth_sL_ip(IJKN,M,L),J)
489                    Kth_sL_s = AVG_Y_S(Kth_sL_ip(IJKS,M,L),Kth_sL_ip(IJK,M,L), JM)
490                    Kth_sL_t = AVG_Z_S(Kth_sL_ip(IJK,M,L), Kth_sL_ip(IJKT,M,L),K)
491                    Kth_sL_b = AVG_Z_S(Kth_sL_ip(IJKB,M,L),Kth_sL_ip(IJK,M,L), KM)
492     
493                    S21a_sum = S21a_sum + ( (Kth_sL_e*(T_PL_E-T_PL_p) )*&
494                         ODX_E(I)*AYZ(IJK) - (Kth_sL_w*(T_PL_p-T_PL_W) )*ODX_E(IM)*&
495                         AYZ(IMJK) )
496                    S21b_sum = S21b_sum + ( (Kth_sL_n*(T_PL_N-T_PL_p) )*&
497                         ODY_N(J)*AXZ(IJK) - (Kth_sL_s*(T_PL_p-T_PL_S) )*ODY_N(JM)*&
498                         AXZ(IJMK) )
499                    S21c_sum = S21c_sum + ( (Kth_sL_t*(T_PL_T-T_PL_p) )*&
500                         ODZ_T(K)*OX(I)*AXY(IJK) - (Kth_sL_b*(T_PL_p-T_PL_B) )*&
501                         ODZ_T(KM)*OX(I)*AXY(IJKM) )
502     
503     ! Part of Heat Flux: div (q)
504     ! Knu_s_ip*[ni*grad(np)-np*grad(ni)]
505     ! Note S14 terms should evaluate to zero for particles from the same phase
506                    Knu_sL_e = AVG_X_S(Knu_sL_ip(IJK,M,L), Knu_sL_ip(IJKE,M,L),I)
507                    Knu_sM_e = AVG_X_S(Knu_sM_ip(IJK,M,L), Knu_sM_ip(IJKE,M,L),I)
508                    Knu_sL_w = AVG_X_S(Knu_sL_ip(IJKW,M,L),Knu_sL_ip(IJK,M,L), IM)
509                    Knu_sM_w = AVG_X_S(Knu_sM_ip(IJKW,M,L),Knu_sM_ip(IJK,M,L), IM)
510                    Knu_sL_n = AVG_Y_S(Knu_sL_ip(IJK,M,L), Knu_sL_ip(IJKN,M,L),J)
511                    Knu_sM_n = AVG_Y_S(Knu_sM_ip(IJK,M,L), Knu_sM_ip(IJKN,M,L),J)
512                    Knu_sL_s = AVG_Y_S(Knu_sL_ip(IJKS,M,L),Knu_sL_ip(IJK,M,L), JM)
513                    Knu_sM_s = AVG_Y_S(Knu_sM_ip(IJKS,M,L),Knu_sM_ip(IJK,M,L), JM)
514                    Knu_sL_t = AVG_Z_S(Knu_sL_ip(IJK,M,L), Knu_sL_ip(IJKT,M,L),K)
515                    Knu_sM_t = AVG_Z_S(Knu_sM_ip(IJK,M,L), Knu_sM_ip(IJKT,M,L),K)
516                    Knu_sL_b = AVG_Z_S(Knu_sL_ip(IJKB,M,L),Knu_sL_ip(IJK,M,L), KM)
517                    Knu_sM_b = AVG_Z_S(Knu_sM_ip(IJKB,M,L),Knu_sM_ip(IJK,M,L), KM)
518     
519                    S14a_sum = S14a_sum + ( (Knu_sL_e*(NU_PL_E-NU_PL_p) - &
520                         Knu_sM_e*(NU_PM_E-NU_PM_p) )*ODX_E(I)*AYZ(IJK) - (Knu_sL_w*&
521                         (NU_PL_p-NU_PL_W) - Knu_sM_w*(NU_PM_p-NU_PM_W) )*ODX_E(IM)*&
522                         AYZ(IMJK) )
523                    S14b_sum = S14b_sum + ( (Knu_sL_n*(NU_PL_N-NU_PL_p) - &
524                         Knu_sM_n*(NU_PM_N-NU_PM_p) )*ODY_N(J)*AXZ(IJK) - (Knu_sL_s*&
525                         (NU_PL_p-NU_PL_S) - Knu_sM_s*(NU_PM_p-NU_PM_S) )*ODY_N(JM)*&
526                         AXZ(IJMK) )
527                    S14c_sum = S14c_sum + ( (Knu_sL_t*(NU_PL_T-NU_PL_p) - &
528                         Knu_sM_t*(NU_PM_T-NU_PM_p) )*ODZ_T(K)*OX(I)*AXY(IJK) - &
529                         (Knu_sL_b*(NU_PL_p-NU_PL_B) - Knu_sM_b*(NU_PM_p-NU_PM_B) )*&
530                         ODZ_T(KM)*OX(I)*AXY(IJKM) )
531     
532     ! Part of Heat Flux: div (q)
533     ! Kvel_s_ip*[vi-vp]
534     ! Note S9 terms should evaluate to zero for particles from the same phase
535                    Kvel_s_e = AVG_X_H(Kvel_s_ip(IJK,M,L), Kvel_s_ip(IJKE,M,L),I)
536                    Kvel_s_w = AVG_X_H(Kvel_s_ip(IJKW,M,L),Kvel_s_ip(IJK,M,L), IM)
537                    Kvel_s_n = AVG_Y_H(Kvel_s_ip(IJK,M,L), Kvel_s_ip(IJKN,M,L),J)
538                    Kvel_s_s = AVG_Y_H(Kvel_s_ip(IJKS,M,L),Kvel_s_ip(IJK,M,L), JM)
539                    Kvel_s_t = AVG_Z_H(Kvel_s_ip(IJK,M,L), Kvel_s_ip(IJKT,M,L),K)
540                    Kvel_s_b = AVG_Z_H(Kvel_s_ip(IJKB,M,L),Kvel_s_ip(IJK,M,L), KM)
541     
542                    S9_sum = S9_sum + ( Kvel_s_e*(UsM_e-UsL_e)*AYZ(IJK) - &
543                         Kvel_s_w*(UsM_w-UsL_w)*AYZ(IMJK) + Kvel_s_n*(VsM_n-VsL_n)*AXZ(IJK)-&
544                         Kvel_s_s*(VsM_s-VsL_s)*AXZ(IJMK) + Kvel_s_t*(WsM_t-WsL_t)*AXY(IJK)-&
545                         Kvel_s_b*(WsM_b-WsL_b)*AXY(IJKM) )
546     
547               ENDIF    ! (IF M.NE.L)
548     
549           ENDDO
550     
551     !  WARNING: The terms due to granular temperature gradients S21 (a,b,c) have caused
552     !           some converegence issues, remove them from LHS and RHS for debugging (sof).
553     
554           SOURCELHS = ( (S11_sum_lhs+S12_sum_lhs)+&
555               (S10_lhs+S16_lhs+S17_sum_lhs+&
556               S18_sum_lhs-S20_sum_lhs+S13_sum_lhs)*VOL(IJK) + &
557               ZMAX(S21a_sum+S21b_sum+S21c_sum)+ &
558               ZMAX(S14a_sum+S14b_sum+S14c_sum)+ ZMAX(S9_sum) ) / &
559               Theta_m(IJK,M)
560     
561           SOURCERHS = ( S10_rhs+S15_rhs+S16_rhs+S17_sum_rhs+S18_sum_rhs+&
562               S13_sum_rhs) * VOL(IJK) + S11_sum_rhs+S12_sum_rhs+ &
563               ZMAX(- (S14a_sum+S14b_sum+S14c_sum) ) + ZMAX(-S9_sum) + &
564               ZMAX(- (S21a_sum+S21b_sum+S21c_sum) )
565     
566     
567           RETURN
568           END SUBROUTINE SOURCE_GRANULAR_ENERGY_IA
569     
570     
571     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
572     !                                                                      C
573     !  Subroutine: SOURCE_GRANULAR_ENERGY_GD                               C
574     !  Purpose: Calculate the source terms of the Mth pahse granular       C
575     !     energy equation                                                  C
576     !                                                                      C
577     !  Author: Janine E. Galvin                                            C
578     !                                                                      C
579     !  Literature/Document References:                                     C
580     !    Garzo, V., and Dufty, J., Homogeneous cooling state for a         C
581     !      granular mixture, Physical Review E, 1999, Vol 60 (5), 5706-    C
582     !      5713                                                            C
583     !    Garzo, Tenneti, Subramaniam, Hrenya, J. Fluid Mech., 2012, 712,   C
584     !      pp 129-404                                                      C
585     !                                                                      C
586     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
587     
588           SUBROUTINE SOURCE_GRANULAR_ENERGY_GD(SOURCELHS, &
589                         SOURCERHS, IJK, M)
590     
591     !-----------------------------------------------
592     !  Modules
593     !-----------------------------------------------
594           USE compar
595           USE constant
596           USE drag
597           USE fldvar
598           USE fun_avg
599           USE functions
600           USE geometry
601           USE indices
602           USE parallel
603           USE param
604           USE param1
605           USE physprop
606           USE rdf
607           USE residual
608           USE run
609           USE toleranc
610           USE trace
611           USE turb
612           USE visc_g
613           USE visc_s
614           use kintheory
615           IMPLICIT NONE
616     !-----------------------------------------------
617     ! Dummy Arguments
618     !-----------------------------------------------
619     ! Source terms to be kept on rhs (source vector)
620           DOUBLE PRECISION, INTENT(INOUT) :: sourcerhs
621     ! Source terms to be kept on lhs (center coefficient)
622           DOUBLE PRECISION, INTENT(INOUT) :: sourcelhs
623     ! Solid phase index
624           INTEGER, INTENT(IN) :: M
625     ! Indices
626           INTEGER, INTENT(IN) :: IJK
627     !-----------------------------------------------
628     ! Local variables
629     !-----------------------------------------------
630     ! Indices
631           INTEGER :: I, J, K, IM, JM, KM, IMJK, IJMK, IJKM,&
632                      IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
633     ! number densities
634           DOUBLE PRECISION :: NU_PM_E, NU_PM_W, NU_PM_N, NU_PM_S, NU_PM_T,&
635                               NU_PM_B, NU_PM_p
636     ! particle characteristics
637           DOUBLE PRECISION :: M_PM, D_PM
638     ! coefficients in heat flux term
639           DOUBLE PRECISION :: Knu_sM_e, Knu_sM_w, Knu_sM_n, Knu_sM_s,&
640                               Knu_sM_t, Knu_sM_b
641     ! Source terms to be kept on rhs
642           DOUBLE PRECISION :: S8_rhs, S9_rhs, S10_rhs, S11_rhs, S12_rhs,&
643                               S13a_rhs, S13b_rhs, S14a_rhs, S14b_rhs, S15a_rhs,&
644                               S15b_rhs
645     ! Source terms to be kept on lhs
646           DOUBLE PRECISION :: S8_lhs, S10_lhs, S11_lhs, S12_lhs, S13a_lhs,&
647                               S13b_lhs, S14a_lhs, S14b_lhs, S15a_lhs, S15b_lhs
648     ! Source terms from interstitial effects
649           DOUBLE PRECISION :: VSLIP, Kslip, Tslip_rhs, Tslip_lhs, Tvis_lhs
650           DOUBLE PRECISION :: ep_sM
651     ! local var. for GTSH theory
652           DOUBLE PRECISION :: chi, Sgama_lhs, Spsi_rhs
653     !-----------------------------------------------
654     ! Function subroutines
655     !-----------------------------------------------
656           DOUBLE PRECISION, EXTERNAL :: G_gtsh
657     !-----------------------------------------------
658     
659           I = I_OF(IJK)
660           J = J_OF(IJK)
661           K = K_OF(IJK)
662           IM = Im1(I)
663           JM = Jm1(J)
664           KM = Km1(K)
665           IMJK = IM_OF(IJK)
666           IJMK = JM_OF(IJK)
667           IJKM = KM_OF(IJK)
668           IJKE = EAST_OF(IJK)
669           IJKW = WEST_OF(IJK)
670           IJKN = NORTH_OF(IJK)
671           IJKS = SOUTH_OF(IJK)
672           IJKT = TOP_OF(IJK)
673           IJKB = BOTTOM_OF(IJK)
674     
675           S8_lhs = ZERO
676           S10_lhs = ZERO
677           S11_lhs = ZERO
678           S12_lhs = ZERO
679           S13a_lhs = ZERO
680           S13b_lhs = ZERO
681           S14a_lhs = ZERO
682           S14b_lhs = ZERO
683           S15a_lhs = ZERO
684           S15b_lhs = ZERO
685           S8_rhs = ZERO
686           S9_rhs = ZERO
687           S10_rhs = ZERO
688           S11_rhs = ZERO
689           S12_rhs = ZERO
690           S13a_rhs = ZERO
691           S13b_rhs = ZERO
692           S14a_rhs = ZERO
693           S14b_rhs = ZERO
694           S15a_rhs = ZERO
695           S15b_rhs = ZERO
696           Sgama_lhs = ZERO  ! for GTSH theory
697           Spsi_rhs = ZERO   ! for GTSH theory
698     
699           D_PM = D_P(IJK,M)
700           M_PM = (Pi/6.d0)*D_PM**3 * RO_S(IJK,M)
701           EP_SM = EP_s(IJK,M)
702           CHI = G_0(IJK,M,M)
703     
704           NU_PM_p = ROP_S(IJK,M)/M_PM
705           NU_PM_E = ROP_S(IJKE,M)/M_PM
706           NU_PM_W = ROP_S(IJKW,M)/M_PM
707           NU_PM_N = ROP_S(IJKN,M)/M_PM
708           NU_PM_S = ROP_S(IJKS,M)/M_PM
709           NU_PM_T = ROP_S(IJKT,M)/M_PM
710           NU_PM_B = ROP_S(IJKB,M)/M_PM
711     
712     ! Production by shear: (S:grad(v))
713     ! P_s*tr(D)
714           S8_lhs = P_S_C(IJK,M) * ZMAX(TRD_S_C(IJK,M))
715           S8_rhs = P_S_C(IJK,M) * ZMAX(-TRD_S_C(IJK,M))
716     
717     ! Production by shear: (S:grad(v))
718     ! Mu_s*tr(D^2)
719           S9_rhs = 2.d0*Mu_s_c(IJK,M)*TRD_S2(IJK,M)
720     
721     ! Production by shear: (S:grad(v))
722     ! Lambda_s*tr(D)^2
723           S10_lhs = (TRD_S_C(IJK,M)**2)*ZMAX( -LAMBDA_s_C(IJK,M) )
724           S10_rhs = (TRD_S_C(IJK,M)**2)*ZMAX(  LAMBDA_s_C(IJK,M) )
725     
726     ! Energy dissipation by collisions: (3/2)*n*kboltz*T*zeta0
727     ! linearized (3/2)*rop_s*T*zeta0
728           S11_lhs = (3.d0/2.d0)*EDT_s_ip(IJK,M,M)
729           S11_rhs = (1.d0/2.d0)*EDT_s_ip(IJK,M,M)*Theta_m(IJK,M)
730     
731     ! Energy dissipation by collisions: (3/2)*n*kboltz*T*zeta1
732     ! (3/2)*rop_s*T*zeta1
733           S12_lhs = ZMAX( EDvel_sM_ip(IJK,M,M) * TRD_S_C(IJK,M) )
734           S12_rhs = ZMAX( -EDvel_sM_ip(IJK,M,M) * TRD_S_C(IJK,M) )*Theta_m(IJK,M)
735     
736     ! for GTSH theory, the dissipation terms above need to be multiplied
737     ! by 3/2 rop_s(ijk,m)
738     ! GTSH theory has two additional terms as source (Psi) and sink (gama)
739     ! in eq (4.9) of GTSH JFM paper
740           IF(KT_TYPE_ENUM == GTSH_2012) THEN
741             S11_lhs = 1.5d0*ROP_s(IJK,M) * S11_lhs
742             S11_rhs = 1.5d0*ROP_s(IJK,M) * S11_rhs
743             S12_lhs = 1.5d0*ROP_s(IJK,M) * S12_lhs
744             S12_rhs = 1.5d0*ROP_s(IJK,M) * S12_rhs
745             Sgama_lhs = 3d0*NU_PM_p * G_gtsh(EP_SM, chi, IJK, M)
746             Spsi_rhs = 1.5d0*ROP_s(IJK,M) * xsi_gtsh(ijk)
747           ENDIF
748     
749     ! Part of Heat Flux: div (q)
750     ! Knu_s_ip*grad(nu)
751           Knu_sM_e = AVG_X_S(Kphi_s(IJK,M), Kphi_s(IJKE,M),I)
752           Knu_sM_w = AVG_X_S(Kphi_s(IJKW,M),Kphi_s(IJK,M), IM)
753           Knu_sM_n = AVG_Y_S(Kphi_s(IJK,M), Kphi_s(IJKN,M),J)
754           Knu_sM_s = AVG_Y_S(Kphi_s(IJKS,M),Kphi_s(IJK,M), JM)
755           Knu_sM_t = AVG_Z_S(Kphi_s(IJK,M), Kphi_s(IJKT,M),K)
756           Knu_sM_b = AVG_Z_S(Kphi_s(IJKB,M),Kphi_s(IJK,M), KM)
757     
758           S13a_rhs = Knu_sM_e*ZMAX( (NU_PM_E-NU_PM_p) )*ODX_E(I)*AYZ(IJK)
759           S13a_lhs = Knu_sM_e*ZMAX( -(NU_PM_E-NU_PM_p) )*ODX_E(I)*AYZ(IJK)
760     
761           S13b_rhs = Knu_sM_w*ZMAX( -(NU_PM_p-NU_PM_W) )*ODX_E(IM)*AYZ(IMJK)
762           S13b_lhs = Knu_sM_w*ZMAX( (NU_PM_p-NU_PM_W) )*ODX_E(IM)*AYZ(IMJK)
763     
764           S14a_rhs = Knu_sM_n*ZMAX( (NU_PM_N-NU_PM_p) )*ODY_N(J)*AXZ(IJK)
765           S14a_lhs = Knu_sM_n*ZMAX( -(NU_PM_N-NU_PM_p) )*ODY_N(J)*AXZ(IJK)
766     
767           S14b_rhs = Knu_sM_s*ZMAX( -(NU_PM_p-NU_PM_S) )*ODY_N(JM)*AXZ(IJMK)
768           S14b_lhs = Knu_sM_s*ZMAX( (NU_PM_p-NU_PM_S) )*ODY_N(JM)*AXZ(IJMK)
769     
770           S15a_rhs = Knu_sM_t*ZMAX( (NU_PM_T-NU_PM_p) )*ODZ_T(K)*OX(I)*AXY(IJK)
771           S15a_lhs = Knu_sM_t*ZMAX( -(NU_PM_T-NU_PM_p) )*ODZ_T(K)*OX(I)*AXY(IJK)
772     
773           S15b_rhs = Knu_sM_b*ZMAX( -(NU_PM_p-NU_PM_B) )*ODZ_T(KM)*OX(I)*AXY(IJKM)
774           S15b_lhs = Knu_sM_b*ZMAX( (NU_PM_p-NU_PM_B) )*ODZ_T(KM)*OX(I)*AXY(IJKM)
775     
776           SOURCELHS = ( (S8_lhs+S10_lhs)*VOL(IJK) + &
777               S13a_lhs+S13b_lhs+S14a_lhs+S14b_lhs+S15a_lhs+S15b_lhs)/Theta_m(IJK,M)&
778               + (S11_lhs + S12_lhs + Sgama_lhs)*VOL(IJK)
779     
780     
781           SOURCERHS = ( S8_rhs+S9_rhs+S10_rhs+S11_rhs+S12_rhs+Spsi_rhs) * VOL(IJK) + &
782               S13a_rhs+S13b_rhs+S14a_rhs+S14b_rhs+S15a_rhs+S15b_rhs
783     
784     
785     ! this is only done for GD_99, do not add these terms to GTSH
786           IF(SWITCH > ZERO .AND. RO_g0 /= ZERO .AND. &
787              (KT_TYPE_ENUM == GD_1999)) THEN
788     
789              VSLIP = (U_S(IJK,M)-U_G(IJK))**2 + (V_S(IJK,M)-V_G(IJK))**2 +&
790                 (W_S(IJK,M)-W_G(IJK))**2
791              VSLIP = DSQRT(VSLIP)
792     
793     ! production by gas-particle slip: Koch & Sangani (1999)
794              Kslip = SWITCH*81.d0*EP_SM*(MU_G(IJK)*VSLIP)**2.d0 / &
795                 (chi*D_P(IJK,M)**3.D0*RO_S(IJK,M)*DSQRT(PI))
796     
797              Tslip_rhs = 1.5d0*Kslip/( (THETA_M(IJK,M)+SMALL_NUMBER)**0.50)*VOL(IJK)
798              Tslip_lhs = 0.5d0*Kslip/( (THETA_M(IJK,M)+SMALL_NUMBER)**1.50)*VOL(IJK)
799     
800     ! dissipation by viscous damping: Gidaspow (1994)
801              Tvis_lhs = SWITCH*3d0*F_GS(IJK,M)*VOL(IJK)
802     
803              SOURCELHS = SOURCELHS + Tslip_lhs + Tvis_lhs
804              SOURCERHS = SOURCERHS + Tslip_rhs
805           ENDIF
806     
807           RETURN
808           END SUBROUTINE SOURCE_GRANULAR_ENERGY_GD
809     
810     
811