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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: SOURCE_GHD_GRANULAR_ENERGY(sourcelhs,sourcerhs,IJK,IER)
4     !  Purpose: Calculate the source terms in the granular energy equation C
5     !           for GHD theory                                             C
6     !                                                                      C
7     !  Author: S. Benyahia, J. Galvin, C. Hrenya        Date: 22-JAN-09    C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Literature/Document References:                                     C
11     !     C. Hrenya handnotes and Garzo, Hrenya, Dufty papers (PRE, 2007)  C
12     !                                                                      C
13     !  Variables referenced:                                               C
14     !                                                                      C
15     !  Variables modified:                                                 C
16     !                                                                      C
17     !     Local variables: sourcelhs, sourcerhs                            C
18     !                                                                      C
19     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20           SUBROUTINE SOURCE_GHD_GRANULAR_ENERGY(SOURCELHS, SOURCERHS, IJK)
21     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
22     !...Switches: -xf
23     !
24     !     Include param.inc file to specify parameter values
25     !
26     !-----------------------------------------------
27     !   M o d u l e s
28     !-----------------------------------------------
29           USE param
30           USE param1
31           USE parallel
32           USE physprop
33           USE run
34           USE drag
35           USE geometry
36           USE fldvar
37           USE visc_s
38           USE ghdtheory
39           USE trace
40           USE indices
41           USE constant
42           USE toleranc
43           USE compar        !//d
44           USE fun_avg
45           USE functions
46           IMPLICIT NONE
47     !-----------------------------------------------
48     !   G l o b a l   P a r a m e t e r s
49     !-----------------------------------------------
50     !-----------------------------------------------
51     !   D u m m y   A r g u m e n t s
52     !-----------------------------------------------
53     !
54     !                      Indices
55           INTEGER          IJK, I, J, K, IM, JM, KM, IMJK, IJMK, IJKM,&
56                            IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
57     !
58     !                      number densities
59           DOUBLE PRECISION NiE, NiW, NiN, NiS, NiT,&
60                            NiB, Nip, Ntotal
61     !
62     !                      mass of particles & non-zero gran. temp.
63           DOUBLE PRECISION Mi, NonZeroTheta
64     !
65     !                      Thermal mobility-related source terms
66           DOUBLE PRECISION ThermMobilityX, ThermMobilityY, ThermMobilityZ
67     !
68     !                      del.Joi and Fi.Joi terms
69           DOUBLE PRECISION DelDotJoi, FiDotJoi, JoiXC, JoiYC, JoiZC
70           DOUBLE PRECISION FiXC, FiYC, FiZC,ni(smax), SIGMAI(smax)
71     !
72     !                      phase index
73           INTEGER          M, L
74     !
75     !                      Dufour-related source terms
76           DOUBLE PRECISION DufourX, DufourY, DufourZ, DijQTerm, &
77                            DijQTermE, DijQTermW, DijQTermN, DijQTermS, DijQTermT, DijQTermB,&
78                            LijTermW,LijTermE,LijTermN,LijTermS,LijTermT,LijTermB
79     
80           DOUBLE PRECISION DijQTermE_H,DijQTermE_A,DijQTermW_H,DijQTermW_A,DijQTermN_H,DijQTermN_A,DijQTermS_H,&
81                            DijQTermS_A,DijQTermT_H,DijQTermT_A,DijQTermB_H,DijQTermB_A,LijTermE_H,LijTermE_A, &
82                            LijTermW_H,LijTermW_A,LijTermN_H,LijTermN_A,LijTermS_H,LijTermS_A,LijTermT_H, &
83                            LijTermT_A,LijTermB_H,LijTermB_A
84     !
85     !                      Source terms to be kept on RHS
86           DOUBLE PRECISION SOURCERHS, PressureRhs, ShearProduction, BulkViscRhs, DissDivURhs, phi_tot, &
87                            SOURCE_FLUID,SINK_FLUID
88     !    DOUBLE PRECISION chi_ij
89     !
90     !                      Source terms to be kept on LHS
91           DOUBLE PRECISION SOURCELHS, PressureLhs, CollDissipation, BulkViscLhs, DissDivULhs,VSLIP
92           DOUBLE PRECISION UGC, USCM, VGC, VSCM, WGC, WSCM
93     
94           I = I_OF(IJK)
95           J = J_OF(IJK)
96           K = K_OF(IJK)
97           IM = Im1(I)
98           JM = Jm1(J)
99           KM = Km1(K)
100           IMJK = IM_OF(IJK)
101           IJMK = JM_OF(IJK)
102           IJKM = KM_OF(IJK)
103           IJKE = EAST_OF(IJK)
104           IJKW = WEST_OF(IJK)
105           IJKN = NORTH_OF(IJK)
106           IJKS = SOUTH_OF(IJK)
107           IJKT = TOP_OF(IJK)
108           IJKB = BOTTOM_OF(IJK)
109     
110           NonZeroTheta = MAX(THETA_M(IJK,MMAX), SMALL_NUMBER)
111     
112           Ntotal = ZERO
113           phi_tot = ZERO
114           DO M = 1,SMAX
115               Ntotal = Ntotal + ROP_S(IJK,M)*6.d0/(Pi*D_P(IJK,M)**3 * RO_S(IJK,M))
116               ni(M) = ROP_S(IJK,M)*6.d0/(Pi*D_P(IJK,M)**3 * RO_S(IJK,M))
117               phi_tot = phi_tot + ROP_S(IJK,M)/RO_S(IJK,M)
118               SIGMAI(M) = D_P(IJK,M)
119           ENDDO
120     
121     ! Production by shear: (S:grad(v))
122     !     p_s*tr(D)
123           PressureRhs = P_S_C(IJK,MMAX)*ZMAX(( -TRD_S_C(IJK,MMAX) ))
124           PressureLhs = P_S_C(IJK,MMAX)*ZMAX((  TRD_S_C(IJK,MMAX) ))
125     
126     !     mu_s*tr(D^2)
127           ShearProduction = 2d0*MU_S_C(IJK,MMAX)*TRD_S2(IJK,MMAX)
128     
129     !     lambda_s*tr(D)^2
130           BulkViscRhs = ZMAX(  LAMBDA_S_C(IJK,MMAX) ) * TRD_S_C(IJK,MMAX)**2
131           BulkViscLhs = ZMAX( -LAMBDA_S_C(IJK,MMAX) ) * TRD_S_C(IJK,MMAX)**2
132     
133     
134     ! Energy dissipation by collisions: (3/2)*n*T*zeta
135     !     (3/2)*n*T*zeta0; zeroth order cooling rate term
136           CollDissipation = 1.5d0*Ntotal*Zeta0(IJK)
137     !     (3/2)*n*T*zetaU*div(U) :
138           DissDivURhs = 1.5d0*Ntotal*Theta_m(IJK,MMAX)* ZMAX( -ZetaU(IJK)*TRD_S_C(IJK,MMAX) )
139           DissDivULhs = 1.5d0*Ntotal* ZMAX(  ZetaU(IJK)*TRD_S_C(IJK,MMAX) )
140     
141     
142           DufourX = ZERO
143           DufourY = ZERO
144           DufourZ = ZERO
145           ThermMobilityX = ZERO
146           ThermMobilityY = ZERO
147           ThermMobilityZ = ZERO
148           DelDotJoi = ZERO
149           FiDotJoi  = ZERO
150           SOURCE_FLUID = ZERO
151           SINK_FLUID = ZERO
152           DO M = 1,SMAX
153     
154     ! Part of heat flux: div (q)
155     !     Sum_ij [ div( T^2*DijQ/nj*grad(nj)) ] -> Dufour term
156     
157               DO L = 1,SMAX
158                   Mi = (Pi/6.d0)*D_P(IJK,L)**3 * RO_S(IJK,L)
159     
160                   Nip = ROP_S(IJK,L) /Mi
161                   NiE = ROP_S(IJKE,L)/Mi
162                   NiW = ROP_S(IJKW,L)/Mi
163                   NiN = ROP_S(IJKN,L)/Mi
164                   NiS = ROP_S(IJKS,L)/Mi
165     
166                   DijQTerm  = zero
167                   DijQTermE = zero
168                   DijQTermW = zero
169                   DijQTermN = zero
170                   DijQTermS = zero
171     
172                   if(ROP_S(IJK,L)/RO_S(IJK,L) > DIL_EP_S) DijQTerm = Theta_m(IJK,MMAX)**2*DijQ(IJK,M,L) / Nip
173                   if(ROP_S(IJKE,L)/RO_S(IJKE,L) > DIL_EP_S) DijQTermE =Theta_m(IJKE,MMAX)**2*DijQ(IJKE,M,L) / NiE
174                   if(ROP_S(IJKW,L)/RO_S(IJKW,L) > DIL_EP_S) DijQTermW =Theta_m(IJKW,MMAX)**2*DijQ(IJKW,M,L) / NiW
175                   if(ROP_S(IJKN,L)/RO_S(IJKN,L) > DIL_EP_S) DijQTermN =Theta_m(IJKN,MMAX)**2*DijQ(IJKN,M,L) / NiN
176                   if(ROP_S(IJKS,L)/RO_S(IJKS,L) > DIL_EP_S) DijQTermS =Theta_m(IJKS,MMAX)**2*DijQ(IJKS,M,L) / NiS
177     
178                     DijQTermE_H = AVG_X_S(DijQTerm, DijQTermE, I)
179                     DijQTermE_A = AVG_X(DijQTerm, DijQTermE, I)
180                   IF(MIN(ABS(DijQTermE_H),ABS(DijQTermE_A)) .eq. ABS(DijQTermE_H))THEN
181                     DijQTermE = DijQTermE_H
182                   ELSE
183                     DijQTermE = DijQTermE_A
184                   ENDIF
185     
186                     DijQTermW_H = AVG_X_S(DijQTermW, DijQTerm, IM)
187                     DijQTermW_A = AVG_X(DijQTermW, DijQTerm, IM)
188                   IF(MIN(ABS(DijQTermW_H),ABS(DijQTermW_A)) .eq. ABS(DijQTermW_H))THEN
189                     DijQTermW = DijQTermW_H
190                   ELSE
191                     DijQTermW = DijQTermW_A
192                   ENDIF
193     
194                     DijQTermN_H = AVG_Y_S(DijQTerm, DijQTermN, J)
195                     DijQTermN_A = AVG_Y(DijQTerm, DijQTermN, J)
196                   IF(MIN(ABS(DijQTermN_H),ABS(DijQTermN_A)) .eq. ABS(DijQTermN_H))THEN
197                     DijQTermN = DijQTermN_H
198                   ELSE
199                     DijQTermN = DijQTermN_A
200                   ENDIF
201     
202                     DijQTermS_H = AVG_Y_S(DijQTermS, DijQTerm, JM)
203                     DijQTermS_A = AVG_Y(DijQTermS, DijQTerm, JM)
204                   IF(MIN(ABS(DijQTermS_H),ABS(DijQTermS_A)) .eq. ABS(DijQTermS_H))THEN
205                     DijQTermS = DijQTermS_H
206                   ELSE
207                     DijQTermS = DijQTermS_A
208                   ENDIF
209     
210                   DufourX = DufourX + ( DijQTermE*(NiE-Nip)*ODX_E(I) - &
211                                         DijQTermW*(Nip-NiW)*ODX_E(IM) )*AYZ(IJK)
212                   DufourY = DufourY + ( DijQTermN*(NiN-Nip)*ODY_N(J) - &
213                                         DijQTermS*(Nip-NiS)*ODY_N(JM) )*AXZ(IJK)
214     
215                   IF(.NOT. NO_K) THEN
216                      NiT = ROP_S(IJKT,L)/Mi
217                      NiB = ROP_S(IJKB,L)/Mi
218     
219                      DijQTermT = zero
220                      DijQTermB = zero
221     
222                      if(ROP_S(IJKT,L)/RO_S(IJKT,L) > DIL_EP_S) DijQTermT = Theta_m(IJK,MMAX)**2*DijQ(IJK,M,L) / NiT
223                      if(ROP_S(IJKB,L)/RO_S(IJKB,L) > DIL_EP_S) DijQTermB = Theta_m(IJK,MMAX)**2*DijQ(IJK,M,L) / NiB
224     
225                      DijQTermT_H = AVG_Z_S(DijQTerm , DijQTermT, K)
226                      DijQTermT_A = AVG_Z(DijQTerm , DijQTermT, K)
227     
228                      IF(MIN(ABS(DijQTermT_H),ABS(DijQTermT_A)) .eq. ABS(DijQTermT_H))THEN
229                        DijQTermT=DijQTermT_H
230                      ELSE
231                        DijQTermT=DijQTermT_A
232                      ENDIF
233     
234                      DijQTermB_H = AVG_Z_S(DijQTermB, DijQTerm, KM)
235                      DijQTermB_A = AVG_Z(DijQTermB, DijQTerm, KM)
236     
237                      IF(MIN(ABS(DijQTermB_H),ABS(DijQTermB_A)) .eq. ABS(DijQTermB_H))THEN
238                        DijQTermB=DijQTermB_H
239                      ELSE
240                        DijQTermB=DijQTermB_A
241                      ENDIF
242     
243                      DufourZ = DufourZ + ( DijQTermT*(NiT-Nip)*ODZ_T(K)*OX(I) - &
244                                            DijQTermB*(Nip-NiB)*ODZ_T(KM)*OX(I) )*AXY(IJK)
245                   ENDIF
246     
247     !     Sum_ij [ div( Lij*Fj) ]; thermal mobility term
248     !     Where Fj = Body Force
249     
250                     LijTermW_H = AVG_X_S(Lij(IJKW,M,L),Lij(IJK,M,L),IM)
251                     LijTermW_A = AVG_X(Lij(IJKW,M,L),Lij(IJK,M,L),IM)
252                   IF(MIN(ABS(LijTermW_H),ABS(LijTermW_A)) .eq. ABS(LijTermW_H))THEN
253                     LijTermW = LijTermW_H
254                   ELSE
255                     LijTermW = LijTermW_A
256                   ENDIF
257     
258                     LijTermE_H = AVG_X_S(Lij(IJK,M,L),Lij(IJKE,M,L),I)
259                     LijTermE_A = AVG_X(Lij(IJK,M,L),Lij(IJKE,M,L),I)
260     
261                   IF(MIN(ABS(LijTermE_H),ABS(LijTermE_A)) .eq. ABS(LijTermE_H))THEN
262                     LijTermE = LijTermE_H
263                   ELSE
264                     LijTermE = LijTermE_A
265                   ENDIF
266     
267                     LijTermN_H = AVG_Y_S(Lij(IJK,M,L),Lij(IJKN,M,L),J)
268                     LijTermN_A = AVG_Y(Lij(IJK,M,L),Lij(IJKN,M,L),J)
269                   IF(MIN(ABS(LijTermN_H),ABS(LijTermN_A)) .eq. ABS(LijTermN_H))THEN
270                     LijTermN = LijTermN_H
271                   ELSE
272                     LijTermN = LijTermN_A
273                   ENDIF
274     
275                     LijTermS_H = AVG_Y_S(Lij(IJKS,M,L),Lij(IJK,M,L),JM)
276                     LijTermS_A = AVG_Y(Lij(IJKS,M,L),Lij(IJK,M,L),JM)
277                   IF(MIN(ABS(LijTermS_H),ABS(LijTermS_A)) .eq. ABS(LijTermS_H))THEN
278                     LijTermS = LijTermS_H
279                   ELSE
280                     LijTermS = LijTermS_A
281                   ENDIF
282     
283                   ThermMobilityX = ThermMobilityX + ( &
284                            FiX(IJK,L) *LijTermE - FiX(IMJK,L)*LijTermW )*AYZ(IJK)
285     
286                   ThermMobilityY = ThermMobilityY + ( &
287                            FiY(IJK,L) *LijTermN - FiY(IJMK,L)*LijTermS )*AXZ(IJK)
288     
289                   IF(.NOT. NO_K) THEN
290     
291                   LijTermT_H = AVG_Z_S(Lij(IJK,M,L),Lij(IJKT,M,L),K)
292                   LijTermT_A = AVG_Z(Lij(IJK,M,L),Lij(IJKT,M,L),K)
293                   IF(MIN(ABS(LijTermT_H),ABS(LijTermT_A)) .eq. ABS(LijTermT_H))THEN
294                     LijTermT = LijTermT_H
295                   ELSE
296                     LijTermT = LijTermT_A
297                   ENDIF
298     
299                   LijTermB_H = AVG_Z_S(Lij(IJKB,M,L),Lij(IJK,M,L),KM)
300                   LijTermB_A = AVG_Z(Lij(IJKB,M,L),Lij(IJK,M,L),KM)
301                   IF(MIN(ABS(LijTermB_H),ABS(LijTermB_A)) .eq. ABS(LijTermB_H))THEN
302                     LijTermB = LijTermB_H
303                   ELSE
304                     LijTermB = LijTermB_A
305                   ENDIF
306     
307                   ThermMobilityZ = ThermMobilityZ + ( &
308                            FiZ(IJK,L) *LijTermT  - FiZ(IJKM,L)*LijTermB )*AXY(IJK)
309                   ENDIF
310               ENDDO ! for L = 1, smax
311     
312     ! Additional term arising from subtraction of 3/2*T*continuity
313     !     + (3/2)*T* Sum_i [ div (joi/mi) ]
314     
315               Mi = (Pi/6.d0)*D_P(IJK,M)**3 * RO_S(IJK,M)
316     
317               DelDotJoi = DelDotJoi + 1.5d0*THETA_M(IJK,MMAX)/Mi * ( &
318                           (JoiX(IJK,M) - JoiX(IMJK,M))*AYZ(IJK) + &
319                           (JoiY(IJK,M) - JoiY(IJMK,M))*AXZ(IJK) + &
320                           (JoiZ(IJK,M) - JoiZ(IJKM,M))*AXY(IJK) )
321     
322     
323     ! Species force dot species mass flux
324     !     Sum_i [ Fi dot joi/mi ]
325     
326     !     Calculate species mass flux components at cell center
327               JoiXC = AVG_X_E(JoiX(IMJK,M),JoiX(IJK,M),I)
328               JoiYC = AVG_Y_N(JoiY(IJMK,M),JoiY(IJK,M))
329               JoiZC = AVG_Z_T(JoiZ(IJKM,M),JoiZ(IJK,M))
330     
331     ! external forces evaluated at cell center
332               FiXC = AVG_X_E(FiX(IMJK,M),FiX(IJK,M),I)
333               FiYC = AVG_Y_N(FiY(IJMK,M),FiY(IJK,M))
334               FiZC = AVG_Z_T(FiZ(IJKM,M),FiZ(IJK,M))
335     
336               FiDotJoi  = FiDotJoi  + ( JoiXC*FiXC + JoiYC*FiYC + JoiZC*FiZC ) / Mi
337     
338     
339               IF(SWITCH > ZERO .AND. RO_g0 /= ZERO) THEN ! do nothing for gran. flow
340                 UGC = AVG_X_E(U_G(IMJK),U_G(IJK),I)
341                 VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
342                 WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
343                 USCM = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
344                 VSCM = AVG_Y_N(V_S(IJMK,M),V_S(IJK,M))
345                 WSCM = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M))
346     
347                 VSLIP = (USCM-UGC)**2 + (VSCM-VGC)**2 + (WSCM-WGC)**2
348                 VSLIP = DSQRT(VSLIP)
349     
350     ! Source/Sink due to fluid do not work well with fluid-solids cases that we run
351     ! uncomment the lines of code below to use (W. Holloway and S. Benyahia).
352     !
353     !            call chi_ij_GHD(smax,M,M,SIGMAi,phi_tot,ni,chi_ij)
354     
355     !            SOURCE_FLUID = SOURCE_FLUID + (81D0*EP_S(IJK,M)*(MU_G(IJK)*VSLIP)**2D0/ &
356     !                   (chi_ij*(D_P(IJK,M)**3D0*RO_S(IJK,M)*THETA_M(IJK,M))**0.5D0))*VOL(IJK)
357     
358     !            SINK_FLUID = SINK_FLUID + 3.d0*F_GS(IJK,M)*THETA_M(IJK,M)/Mi
359               ENDIF
360     
361           ENDDO ! for M = 1, smax
362     
363           SINK_FLUID = SINK_FLUID/NonZeroTheta
364     
365           SOURCERHS = (PressureRhs + ShearProduction + BulkViscRhs + DissDivURhs)*VOL(IJK) &
366                      + ZMAX(DufourX)+ZMAX(DufourY)+ZMAX(DufourZ) &
367                      + ZMAX(ThermMobilityX)+ZMAX(ThermMobilityY)+ZMAX(ThermMobilityZ)      &
368                      + ZMAX(DelDotJoi) + ZMAX(FiDotJoi)*VOL(IJK)
369     
370           SOURCERHS = SOURCERHS + SOURCE_FLUID
371     
372           SOURCELHS = ( (PressureLhs + BulkViscLhs)/NonZeroTheta   + &
373                       (CollDissipation + DissDivULhs + SINK_FLUID) ) * VOL(IJK) + &
374                       ( ZMAX(-DufourX)+ZMAX(-DufourY)+ZMAX(-DufourZ) + &
375                        ZMAX(-ThermMobilityX)+ZMAX(-ThermMobilityY)+ZMAX(-ThermMobilityZ) + &
376                        ZMAX(-DelDotJoi) + ZMAX(-FiDotJoi)*VOL(IJK) )/ NonZeroTheta
377     
378           RETURN
379           END SUBROUTINE SOURCE_GHD_GRANULAR_ENERGY
380