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