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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: DES_GRANULAR_TEMPERATURE                                C
4     !  Purpose: Calculate the DES granular temperature                     C
5     !                                                                      C
6     !                                                                      C
7     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  C
8     !  Reviewer:                                          Date:            C
9     !  Revision: For parallel processing modifications are made to         C
10     !            accomodate ghost cells  (Pradeep G.)                      C
11     !                                                                      C
12     !                                                                      C
13     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
14     
15           SUBROUTINE DES_GRANULAR_TEMPERATURE
16     
17     !-----------------------------------------------
18     ! Modules
19     !-----------------------------------------------
20           USE discretelement
21           USE param
22           USE param1
23           USE run
24           USE geometry
25           USE indices
26           USE compar
27           USE sendrecv
28           USE physprop
29           USE fun_avg
30           USE functions
31           USE mpi_utility
32           IMPLICIT NONE
33     !-----------------------------------------------
34     ! Local Variables
35     !-----------------------------------------------
36     ! indices
37           INTEGER :: I, J, K, IJK
38     !
39           INTEGER :: M, LL
40     ! counter for no. of particles in phase m in cell ijk
41           INTEGER :: NP_PHASE(DIMENSION_3, DES_MMAX)
42     ! temporary variable for mth phase granular temperature in cell ijk
43           DOUBLE PRECISION :: TEMP(DIMENSION_3, DES_MMAX)
44     ! accounted for particles
45           INTEGER :: PC
46     ! squared particle velocity v.v
47           DOUBLE PRECISION :: SQR_VEL, SQR_ROT_VEL
48     !-----------------------------------------------
49     
50     ! Calculate a local species granular temperature for current instant of
51     ! time.  Note that the following calculation of species granular
52     ! temperature employs a fluctuating particle velocity that is defined
53     ! as the difference between a particles velocity and the corresponding
54     ! local mean velocity of that particles species evaluated at the same
55     ! instant of time.
56     !----------------------------------------------------------------->>>
57     ! The following calculations are performed on the 'fluid' grid
58           TEMP(:,:) = ZERO
59           NP_PHASE(:,:) = ZERO
60           PC = 0
61           DO LL = 1, MAX_PIP
62     ! skipping particles that do not exist
63              IF(.NOT.PEA(LL,1)) CYCLE
64              PC = PC + 1
65     ! skipping ghost particles
66              IF(PEA(LL,4)) CYCLE
67     
68              I = PIJK(LL,1)
69              J = PIJK(LL,2)
70              K = PIJK(LL,3)
71              IJK = PIJK(LL,4)
72     
73              M = PIJK(LL,5)
74              NP_PHASE(IJK,M) = NP_PHASE(IJK,M) + 1
75     
76              TEMP(IJK,M) = TEMP(IJK,M) + &
77                 (DES_VEL_NEW(1,LL)-DES_U_s(IJK,M))**2
78              TEMP(IJK,M) = TEMP(IJK,M) + &
79                 (DES_VEL_NEW(2,LL)-DES_V_s(IJK,M))**2
80              IF(DO_K) THEN
81                 TEMP(IJK,M) = TEMP(IJK,M) + &
82                    (DES_VEL_NEW(3,LL)-DES_W_s(IJK,M))**2
83              ENDIF
84     
85              IF(PC .EQ. PIP) EXIT
86           ENDDO
87     
88     ! loop over all fluid cells
89           DO IJK = IJKSTART3, IJKEND3
90              IF(FLUID_AT(IJK)) THEN
91     
92                 DO M = 1,DES_MMAX
93                    IF (NP_PHASE(IJK,M) > 0 ) THEN
94                       DES_THETA(IJK,M) = TEMP(IJK,M)/&
95                          DBLE(DIMN*NP_PHASE(IJK,M))
96                    ELSE
97                       DES_THETA(IJK,M) = ZERO
98                    ENDIF
99                 ENDDO
100              ENDIF
101           ENDDO
102     
103     
104     ! Calculate global quantities: granular temperature, kinetic energy,
105     ! potential energy and average velocity at the current instant of time
106     !----------------------------------------------------------------->>>
107     
108     ! initialization for calculations
109           DES_KE = ZERO
110           DES_ROTE = ZERO
111           DES_PE = ZERO
112           DES_VEL_AVG(:) = ZERO
113     
114     ! Calculate global average velocity, kinetic energy &
115     ! potential energy
116           PC = 0
117           DO LL = 1, MAX_PIP
118              IF(PEA(LL,1)) PC = PC + 1
119     ! skipping ghost particles and particles that don't exist
120              IF(.NOT.PEA(LL,1) .OR. PEA(LL,4)) CYCLE
121     
122              SQR_VEL = ZERO
123              SQR_ROT_VEL = ZERO
124              DO I = 1, DIMN
125                 SQR_VEL = SQR_VEL + DES_VEL_NEW(I,LL)**2
126              ENDDO
127     
128              DO I = 1, merge(1,3,NO_K)
129                 SQR_ROT_VEL = SQR_ROT_VEL + OMEGA_NEW(I,LL)**2
130              ENDDO
131     
132     
133              DES_KE = DES_KE + PMASS(LL)/2.d0 * SQR_VEL
134     ! Calculation of rotational kinetic energy
135              DES_ROTE = DES_ROTE +                                         &
136                 (0.4D0*PMASS(LL)*DES_RADIUS(LL)**2)/2.d0 * SQR_ROT_VEL
137              DES_PE = DES_PE + PMASS(LL)*DBLE(ABS(GRAV(2)))*&
138                 DES_POS_NEW(2,LL)
139              DES_VEL_AVG(:) =  DES_VEL_AVG(:) + DES_VEL_NEW(:,LL)
140     
141              IF(PC .EQ. PIP) EXIT
142           ENDDO
143     
144     !Calculating total number of particles in the entire domain
145     !Needed for correct average velocities and granular temp.
146           CALL GLOBAL_ALL_SUM(PIP-iGHOST_CNT,TOT_PAR)
147           CALL GLOBAL_ALL_SUM(DES_VEL_AVG(1:DIMN))
148     
149     !J.Musser changed PARTICLES TO PIP
150           IF(TOT_PAR > 0) DES_VEL_AVG(:) = DES_VEL_AVG(:)/DBLE(TOT_PAR)
151     
152     ! The following quantities are primarily used for debugging/developing
153     ! and allow a quick check of the energy conservation in the system.
154     ! In their current form they are best applied to monodisperse cases.
155     ! Calculate x,y,z components of global energy & granular temperature
156           GLOBAL_GRAN_ENERGY(:) = ZERO
157           GLOBAL_GRAN_TEMP(:)  = ZERO
158           PC = 0
159           DO LL = 1, MAX_PIP
160              IF(PEA(LL,1)) PC = PC + 1
161     ! skipping ghost particles and particles that don't exist
162              IF(.NOT.PEA(LL,1) .OR. PEA(LL,4)) CYCLE
163     
164              GLOBAL_GRAN_ENERGY(:) = GLOBAL_GRAN_ENERGY(:) + &
165                 0.5d0*PMASS(LL)*(DES_VEL_NEW(:,LL)-DES_VEL_AVG(:))**2
166              GLOBAL_GRAN_TEMP(:) = GLOBAL_GRAN_TEMP(:) + &
167                 (DES_VEL_NEW(:,LL)-DES_VEL_AVG(:))**2
168     
169              IF(PC .EQ. PIP) EXIT
170           ENDDO
171     
172           CALL GLOBAL_ALL_SUM(GLOBAL_GRAN_TEMP)
173           CALL GLOBAL_ALL_SUM(GLOBAL_GRAN_ENERGY)
174           CALL GLOBAL_ALL_SUM(DES_KE)
175           CALL GLOBAL_ALL_SUM(DES_PE)
176           CALL GLOBAL_ALL_SUM(DES_ROTE)
177     
178           IF(TOT_PAR > 0) GLOBAL_GRAN_ENERGY(:) =                          &
179              GLOBAL_GRAN_ENERGY(:)/DBLE(TOT_PAR)
180           IF(TOT_PAR > 0) GLOBAL_GRAN_TEMP(:) =                            &
181              GLOBAL_GRAN_TEMP(:)/DBLE(TOT_PAR)
182     
183     
184     
185           RETURN
186     
187           END SUBROUTINE DES_GRANULAR_TEMPERATURE
188     
189     
190     
191     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
192     !                                                                      !
193     !  Subroutine: CALC_DES_BEDHEIGHT                                      !
194     !  Purpose: Calculate an average bed height for each solids phase      !
195     !           present                                                    !
196     !                                                                      !
197     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
198     
199           SUBROUTINE CALC_DES_BEDHEIGHT
200     
201     !-----------------------------------------------
202     ! Modules
203     !-----------------------------------------------
204           USE indices
205           USE geometry
206           USE compar
207           USE discretelement
208           USE des_bc
209           USE physprop
210           USE fldvar
211           USE functions
212           IMPLICIT NONE
213     !-----------------------------------------------
214     ! Local Variables
215     !-----------------------------------------------
216     ! particle index
217           INTEGER :: L
218     ! accounted for particles
219           INTEGER :: PC
220     ! solids phase no.
221           INTEGER :: M
222     ! ijk indices
223           INTEGER :: J, IJK
224     ! average height of fluid cell
225           DOUBLE PRECISION :: hcell
226     ! height of particle (y-position)
227           DOUBLE PRECISION :: hpart
228     ! volume fraction of phase M in fluid cell
229           DOUBLE PRECISION :: EP_SM
230     ! tmp variables for calculations
231           DOUBLE PRECISION :: tmp_num(DES_MMAX), tmp_den(DES_MMAX)
232     !-----------------------------------------------
233     
234     ! calculation of bed height following the formulation of Goldschmidt et
235     ! al., Powder Technology 138, 2003, 135-159
236           tmp_num(:) = ZERO
237           tmp_den(:) = ZERO
238           DO IJK = IJKSTART3, IJKEND3
239              J = J_OF(IJK)
240              DO M = 1, DES_MMAX
241                 IF(DES_ROP_S(IJK,M) > ZERO) THEN
242                    hcell = 0.5d0*(YN(J)+YN(J-1))
243                    EP_SM = DES_ROP_S(IJK,M)/DES_RO_S(M)
244                    tmp_num(M) = tmp_num(M) + EP_SM*hcell*VOL(IJK)
245                    tmp_den(M) = tmp_den(M) + EP_SM*VOL(IJK)
246                 ENDIF
247              ENDDO
248           ENDDO
249     ! calculate avg height for each phase
250           IF (PIP >0) bed_height(:) = tmp_num(:)/tmp_den(:)
251     
252     
253     ! alternative method to calculating bed height (turned off atm)
254           IF(.FALSE.) THEN
255           tmp_num(:) = ZERO
256           tmp_den(:) = ZERO
257     
258           PC = 0
259           DO L = 1, MAX_PIP
260              IF(PEA(L,1)) PC = PC + 1
261     ! skipping ghost particles and particles that don't exist
262              IF(.NOT.PEA(L,1) .OR. PEA(L,4)) CYCLE
263              M = PIJK(L,5)
264              hpart = DES_POS_NEW(2,L)
265              tmp_num(M) = tmp_num(M) + hpart
266              tmp_den(M) = tmp_den(M) + 1
267              IF(PC .EQ. PIP) EXIT
268           ENDDO
269           ! calculate avg height for each phase
270           IF (PIP >0) bed_height(:) = tmp_num(:)/tmp_den(:)
271           ENDIF
272     
273           RETURN
274     
275           END SUBROUTINE CALC_DES_BEDHEIGHT
276     
277