File: RELATIVE:/../../../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(:,:) = 0
60           PC = 0
61           DO LL = 1, MAX_PIP
62     ! skipping particles that do not exist
63              IF(IS_NONEXISTENT(LL)) CYCLE
64              PC = PC + 1
65     ! skipping ghost particles
66              IF(IS_GHOST(LL) .or. IS_ENTERING_GHOST(LL) .or. IS_EXITING_GHOST(LL)) 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     ! skipping ghost particles and particles that don't exist
119              IF(IS_NONEXISTENT(LL)) CYCLE
120              PC = PC + 1
121              IF(IS_GHOST(LL) .OR. IS_ENTERING_GHOST(LL) .OR. IS_EXITING_GHOST(LL)) CYCLE
122     
123              SQR_VEL = ZERO
124              SQR_ROT_VEL = ZERO
125              DO I = 1, DIMN
126                 SQR_VEL = SQR_VEL + DES_VEL_NEW(I,LL)**2
127              ENDDO
128     
129              DO I = 1, merge(1,3,NO_K)
130                 SQR_ROT_VEL = SQR_ROT_VEL + OMEGA_NEW(I,LL)**2
131              ENDDO
132     
133     
134              DES_KE = DES_KE + PMASS(LL)/2.d0 * SQR_VEL
135     ! Calculation of rotational kinetic energy
136              DES_ROTE = DES_ROTE +                                         &
137                 (0.4D0*PMASS(LL)*DES_RADIUS(LL)**2)/2.d0 * SQR_ROT_VEL
138              DES_PE = DES_PE + PMASS(LL)*DBLE(ABS(GRAV(2)))*&
139                 DES_POS_NEW(2,LL)
140              DES_VEL_AVG(:) =  DES_VEL_AVG(:) + DES_VEL_NEW(:,LL)
141     
142              IF(PC .EQ. PIP) EXIT
143           ENDDO
144     
145     !Calculating total number of particles in the entire domain
146     !Needed for correct average velocities and granular temp.
147           CALL GLOBAL_ALL_SUM(PIP-iGHOST_CNT,TOT_PAR)
148           CALL GLOBAL_ALL_SUM(DES_VEL_AVG(1:DIMN))
149     
150     !J.Musser changed PARTICLES TO PIP
151           IF(TOT_PAR > 0) DES_VEL_AVG(:) = DES_VEL_AVG(:)/DBLE(TOT_PAR)
152     
153     ! The following quantities are primarily used for debugging/developing
154     ! and allow a quick check of the energy conservation in the system.
155     ! In their current form they are best applied to monodisperse cases.
156     ! Calculate x,y,z components of global energy & granular temperature
157           GLOBAL_GRAN_ENERGY(:) = ZERO
158           GLOBAL_GRAN_TEMP(:)  = ZERO
159           PC = 0
160           DO LL = 1, MAX_PIP
161     
162     ! skipping ghost particles and particles that don't exist
163              IF(IS_NONEXISTENT(LL)) CYCLE
164              PC = PC + 1
165              IF(IS_GHOST(LL) .OR. IS_ENTERING_GHOST(LL) .OR. IS_EXITING_GHOST(LL)) CYCLE
166     
167              GLOBAL_GRAN_ENERGY(:) = GLOBAL_GRAN_ENERGY(:) + &
168                 0.5d0*PMASS(LL)*(DES_VEL_NEW(:,LL)-DES_VEL_AVG(:))**2
169              GLOBAL_GRAN_TEMP(:) = GLOBAL_GRAN_TEMP(:) + &
170                 (DES_VEL_NEW(:,LL)-DES_VEL_AVG(:))**2
171     
172              IF(PC .EQ. PIP) EXIT
173           ENDDO
174     
175           CALL GLOBAL_ALL_SUM(GLOBAL_GRAN_TEMP)
176           CALL GLOBAL_ALL_SUM(GLOBAL_GRAN_ENERGY)
177           CALL GLOBAL_ALL_SUM(DES_KE)
178           CALL GLOBAL_ALL_SUM(DES_PE)
179           CALL GLOBAL_ALL_SUM(DES_ROTE)
180     
181           IF(TOT_PAR > 0) GLOBAL_GRAN_ENERGY(:) =                          &
182              GLOBAL_GRAN_ENERGY(:)/DBLE(TOT_PAR)
183           IF(TOT_PAR > 0) GLOBAL_GRAN_TEMP(:) =                            &
184              GLOBAL_GRAN_TEMP(:)/DBLE(TOT_PAR)
185     
186           RETURN
187     
188           END SUBROUTINE DES_GRANULAR_TEMPERATURE
189     
190     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
191     !                                                                      !
192     !  Subroutine: CALC_DES_BEDHEIGHT                                      !
193     !  Purpose: Calculate an average bed height for each solids phase      !
194     !           present                                                    !
195     !                                                                      !
196     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
197     
198           SUBROUTINE CALC_DES_BEDHEIGHT
199     
200     !-----------------------------------------------
201     ! Modules
202     !-----------------------------------------------
203           USE indices
204           USE geometry
205           USE compar
206           USE discretelement
207           USE des_bc
208           USE physprop
209           USE fldvar
210           USE functions
211           IMPLICIT NONE
212     !-----------------------------------------------
213     ! Local Variables
214     !-----------------------------------------------
215     ! particle index
216           INTEGER :: L
217     ! accounted for particles
218           INTEGER :: PC
219     ! solids phase no.
220           INTEGER :: M
221     ! ijk indices
222           INTEGER :: J, IJK
223     ! average height of fluid cell
224           DOUBLE PRECISION :: hcell
225     ! height of particle (y-position)
226           DOUBLE PRECISION :: hpart
227     ! volume fraction of phase M in fluid cell
228           DOUBLE PRECISION :: EP_SM
229     ! tmp variables for calculations
230           DOUBLE PRECISION :: tmp_num(DES_MMAX), tmp_den(DES_MMAX)
231     !-----------------------------------------------
232     
233     ! calculation of bed height following the formulation of Goldschmidt et
234     ! al., Powder Technology 138, 2003, 135-159
235           tmp_num(:) = ZERO
236           tmp_den(:) = ZERO
237           DO IJK = IJKSTART3, IJKEND3
238              J = J_OF(IJK)
239              DO M = 1, DES_MMAX
240                 IF(DES_ROP_S(IJK,M) > ZERO) THEN
241                    hcell = 0.5d0*(YN(J)+YN(J-1))
242                    EP_SM = DES_ROP_S(IJK,M)/DES_RO_S(M)
243                    tmp_num(M) = tmp_num(M) + EP_SM*hcell*VOL(IJK)
244                    tmp_den(M) = tmp_den(M) + EP_SM*VOL(IJK)
245                 ENDIF
246              ENDDO
247           ENDDO
248     ! calculate avg height for each phase
249           IF (PIP >0) bed_height(:) = tmp_num(:)/tmp_den(:)
250     
251     ! alternative method to calculating bed height (turned off atm)
252           IF(.FALSE.) THEN
253           tmp_num(:) = ZERO
254           tmp_den(:) = ZERO
255     
256           PC = 0
257           DO L = 1, MAX_PIP
258     ! skipping ghost particles and particles that don't exist
259              IF(IS_NONEXISTENT(L)) CYCLE
260              PC = PC + 1
261              IF(IS_GHOST(L) .OR. IS_ENTERING_GHOST(L) .OR. IS_EXITING_GHOST(L)) CYCLE
262     
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