File: RELATIVE:/../../../mfix.git/model/calc_p_star.f

1     ! TODO:
2     !   p_star calculation should be based on the sum of volume fractions of
3     !   close-packed solids.
4     
5     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
6     !                                                                      C
7     !  Subroutine: CALC_P_star                                             C
8     !  Purpose: Calculate P_star in cells where solids continuity is       C
9     !     solved                                                           C
10     !                                                                      C
11     !  Author: M. Syamlal                                 Date: 21-AUG-96  C
12     !  Reviewer:                                          Date:            C
13     !                                                                      C
14     !                                                                      C
15     !  Literature/Document References:                                     C
16     !                                                                      C
17     !  Variables referenced:                                               C
18     !  Variables modified: P_STAR,                                         C
19     !     if yu_standish or fedors_landel: ep_star_array,                  C
20     !                                      ep_g_blend_start,               C
21     !                                      ep_g_blend_end                  C
22     !  Local variables:                                                    C
23     !                                                                      C
24     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
25     
26           SUBROUTINE CALC_P_STAR(EP_G, P_STAR)
27     
28     !-----------------------------------------------
29     ! Modules
30     !-----------------------------------------------
31           USE param
32           USE param1
33           USE parallel
34           USE geometry
35           USE indices
36           USE physprop
37           USE constant
38           USE pgcor
39           USE pscor
40           USE ur_facs
41           USE residual
42           USE compar
43           USE run
44           USE visc_s
45           USE solids_pressure
46           USE functions
47           IMPLICIT NONE
48     !-----------------------------------------------
49     ! Dummy arguments
50     !-----------------------------------------------
51     ! Gas volume fraction
52           DOUBLE PRECISION, INTENT(IN) :: EP_g(DIMENSION_3)
53     ! Solids pressure
54           DOUBLE PRECISION, INTENT(INOUT) :: P_star(DIMENSION_3)
55     !-----------------------------------------------
56     ! Local variables
57     !-----------------------------------------------
58     !!   HPF$ align P_star(:) with TT(:)
59     !!   HPF$ align EP_g(:) with TT(:)
60     
61     ! Indices
62           INTEGER :: IJK
63     ! Blend factor
64           DOUBLE PRECISION :: blend
65     !-----------------------------------------------
66     ! External functions
67     !-----------------------------------------------
68           DOUBLE PRECISION :: CALC_EP_STAR
69     !-----------------------------------------------
70     
71     !!$omp parallel do private(ijk)
72     !!   HPF$ independent
73     
74           DO IJK = ijkstart3, ijkend3
75              IF (FLUID_AT(IJK)) THEN
76     
77     
78                 IF (YU_STANDISH .OR. FEDORS_LANDEL) THEN
79     ! if Yu_Standish or Fedors_Landel correlations are used, then
80     ! ep_star_array is modified. this is the only time ep_star_array is
81     ! modified (see set_constprop).  (sof Nov-16-2005)
82                    EP_star_array(ijk) = calc_ep_star(ijk)
83     
84     ! now the values of ep_g_blend_start and ep_g_blend_end need to be
85     ! reassigned based on the new values of ep_star_array
86                    IF(BLENDING_STRESS.AND.TANH_BLEND) THEN
87                       ep_g_blend_start(ijk) = ep_star_array(ijk) * 0.99d0
88                       ep_g_blend_end(ijk)   = ep_star_array(ijk) * 1.01d0
89                    ELSEIF(BLENDING_STRESS.AND.SIGM_BLEND) THEN
90                       ep_g_blend_start(ijk) = EP_star_array(ijk) * 0.97d0
91                       ep_g_blend_end(ijk) = EP_star_array(ijk) * 1.01d0
92                    ELSE
93                       ep_g_blend_start(ijk) = ep_star_array(ijk)
94                       ep_g_blend_end(ijk)   = ep_star_array(ijk)
95                    ENDIF
96                 ENDIF
97     
98                 IF (EP_G(IJK) < EP_g_blend_end(ijk)) THEN
99                    P_STAR(IJK) = NEG_H(EP_G(IJK),EP_g_blend_end(ijk))
100                    IF(BLENDING_STRESS) THEN
101                       blend =  blend_function(IJK)
102                       P_STAR (IJK) = (1.0d0-blend) * P_STAR (IJK)
103                    ENDIF
104                 ELSE
105                    P_STAR(IJK) = ZERO
106                 ENDIF
107              ENDIF
108           ENDDO
109     
110           RETURN
111           END SUBROUTINE CALC_P_STAR
112     
113     
114     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
115     !                                                                      C
116     !  FUNCTION: CALC_ep_star                                              C
117     !  Purpose: calculate the local value of maximum packing               C
118     !                                                                      C
119     !  Author: D. Gera and M. Syamlal                     Date: 31-DEC-02  C
120     !  Reviewer:                                          Date:            C
121     !  Modified: S. Benyahia                              Date: 02-May-05  C
122     !                                                                      C
123     !  Literature/Document References:                                     C
124     !    A.B. Yu and N. Standish. Powder Tech, 52 (1987) 233-241           C
125     !    R.F. Fedors and R.F. Landel. Powder Tech, 23 (1979) 225-231       C
126     !                                                                      C
127     !  Variables referenced:                                               C
128     !  Variables modified:                                                 C
129     !  Local variables:                                                    C
130     !                                                                      C
131     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
132     
133           DOUBLE PRECISION FUNCTION CALC_ep_star(IJK)
134     
135     !-----------------------------------------------
136     ! Modules
137     !-----------------------------------------------
138           USE param
139           USE param1
140           USE fldvar
141           USE geometry
142           USE indices
143           USE physprop
144           USE constant
145           USE toleranc
146           USE compar
147           USE run
148           USE fun_avg
149           USE functions
150           IMPLICIT NONE
151     !-----------------------------------------------
152     ! Dummy arguments
153     !-----------------------------------------------
154     ! IJK index
155           INTEGER, INTENT(IN) :: IJK
156     !-----------------------------------------------
157     ! Local variables
158     !-----------------------------------------------
159     ! Indices
160           INTEGER :: I, J
161     
162     !      DOUBLE PRECISION :: xbar
163     
164     ! start sof modifications (02-May-05)
165     ! maximum packing for the mixture
166            DOUBLE PRECISION :: P_IT(MMAX)
167     ! true maximum packing for the mixture
168            DOUBLE PRECISION :: EPs_max_local
169     ! maximum packing fraction for a binary mixture
170            DOUBLE PRECISION :: P_IJ(MMAX, MMAX)
171     ! particle diameter ratio
172            DOUBLE PRECISION :: R_IJ(MMAX, MMAX)
173     ! fractional solids volume corresponding to P_IJ
174            DOUBLE PRECISION :: X_IJ(MMAX, MMAX)
175     ! fractional solids volume in a mixture
176     ! this is Xj in eq. 22 of Yu-Standish
177            DOUBLE PRECISION :: COMP_X_I(MMAX), SUM_LOCAL
178     ! local aliases for particle diameter, solids volume fraction and the
179     ! maximum solids volume fraction which are used to rearrange solids
180     ! phases from coarsest to finest
181            DOUBLE PRECISION :: DP_TMP(MMAX), EPs_TMP(MMAX), &
182                                EPs_max_TMP(MMAX), old_value
183     !-----------------------------------------------
184     
185     
186           IF (CALL_DQMOM) THEN
187     ! sort particles to start from coarsest to finest particles
188     ! assigning values to local aliases
189              DO I = 1, SMAX
190                 DP_TMP(I) = D_P(IJK,I)
191                 EPs_TMP(I) = EP_s(IJK,I)
192                 EPs_max_TMP(I) = ep_s_max(I)
193              ENDDO
194     ! sorting particles from coarse to fine
195              DO I = 1, SMAX
196                 DO J = I , SMAX
197     ! check if phase J is larger than phase I
198                    IF(DP_TMP(I) < DP_TMP(J)) THEN
199     ! temporarily store phase i diameter
200                       old_value = DP_TMP(I)
201     ! overwrite phase i diameter with smaller phase j diameter
202                       DP_TMP(I) = DP_TMP(J)
203     ! overwrite phase j diameter with stoired phase i diameter
204                       DP_TMP(J) = old_value
205     
206                       old_value = EPs_TMP(I)
207                       EPs_TMP(I) = EPs_TMP(J)
208                       EPs_TMP(J) = old_value
209     
210                       old_value = EPs_max_TMP(I)
211                       EPs_max_TMP(I) = EPs_max_TMP(J)
212                       EPs_max_TMP(J) = old_value
213                    ENDIF
214                 ENDDO
215              ENDDO
216     
217           ELSE  ! not dqmom
218     
219     ! assigning values to local aliases
220              DO I = 1, SMAX
221                 DP_TMP(I) = D_P(IJK,M_MAX(I))
222                 EPs_TMP(I) = EP_s(IJK,M_MAX(I))
223                 EPs_max_TMP(I) = ep_s_max(M_MAX(I))
224              ENDDO
225           ENDIF   ! end if/else (call_dqmom)
226     
227     
228     ! this is the way the algorithm was written by Yu and Standish (sof).
229     ! compute equations 25 in Yu-Standish
230     ! (this is also needed by Fedors_Landel)
231           DO I = 1, SMAX
232              SUM_LOCAL = ZERO
233              DO J = 1, SMAX
234                 IF(I .GE. J) THEN
235                    R_IJ(I,J) = DP_TMP(I)/DP_TMP(J)
236                 ELSE
237                    R_IJ(I,J) = DP_TMP(J)/DP_TMP(I)
238                 ENDIF
239                 SUM_LOCAL = SUM_LOCAL + EPs_TMP(J)
240              ENDDO   ! end do (j=1,smax)
241     
242              IF(SUM_LOCAL > DIL_EP_s) THEN
243     ! fractional solids volume see eq. 20
244                 COMP_X_I(I) = EPs_TMP(I)/SUM_LOCAL
245              ELSE
246     ! return first phase ep_s_max in case very dilute
247                 CALC_EP_star = ONE - EPs_max_TMP(1)
248                 RETURN
249              ENDIF
250           ENDDO   ! end do (i=1,smax)
251     
252     ! Begin YU_STANDISH section
253     ! ---------------------------------------------------------------->>>
254           IF(YU_STANDISH) THEN
255     ! compute equation 23-24 in Yu-Standish
256              DO I = 1, SMAX
257                 DO J = 1, SMAX
258                    IF(R_IJ(I,J) .LE. 0.741d0) THEN
259                       IF(J .LT. I) THEN
260                          X_IJ(I,J) = (ONE - R_IJ(I,J)*R_IJ(I,J))/&
261                                      (2.0d0 -  EPs_max_TMP(I))
262                       ELSE
263                          X_IJ(I,J) = ONE - (ONE - R_IJ(I,J)*R_IJ(I,J))/&
264                                      (2.0d0 -  EPs_max_TMP(I))
265                      ENDIF
266                      P_IJ(I, J) = EPs_max_TMP(I) + EPs_max_TMP(I)*&
267                         (ONE-EPs_max_TMP(I)) * (ONE - 2.35d0*R_IJ(I,J) + &
268                         1.35d0*R_IJ(I,J)*R_IJ(I,J) )
269                    ELSE
270                       P_IJ(I, J) = EPs_max_TMP(I)
271                    ENDIF
272                 ENDDO   ! end do (j=1,smax)
273              ENDDO   ! end do (i=1,smax)
274     
275     ! Compute equation 22
276              EPs_max_local = ONE
277              DO I = 1, SMAX
278                 SUM_LOCAL = ZERO
279     
280                 IF(I .GE. 2) THEN
281                    DO J = 1, (I-1)
282                       IF(P_IJ(I,J) == EPs_max_TMP(I)) THEN
283                          SUM_LOCAL = SUM_LOCAL
284                       ELSE
285                          SUM_LOCAL = SUM_LOCAL + (ONE - EPs_max_TMP(I)/&
286                             P_IJ(I,J))*COMP_X_I(J)/X_IJ(I,J)
287                       ENDIF
288                    ENDDO
289                 ENDIF
290     
291                 IF((I+1) .LE. SMAX) THEN
292                    DO J = (I+1), SMAX
293                       IF( P_IJ(I, J) == EPs_max_TMP(I) ) THEN
294                          SUM_LOCAL = SUM_LOCAL
295                       ELSE
296                          SUM_LOCAL = SUM_LOCAL + (ONE - EPs_max_TMP(I)/&
297                             P_IJ(I, J))*COMP_X_I(J)/X_IJ(I, J)
298                       ENDIF
299                    ENDDO
300                 ENDIF
301     
302                 IF (SUM_LOCAL .NE. ZERO) THEN
303                    P_IT(I) = EPs_max_TMP(I)/(ONE - SUM_LOCAL)
304                 ELSE
305     ! do nothing if particles have same diameter
306                    P_IT(I) = ONE
307                 ENDIF
308     
309                 EPs_max_local = MIN(P_IT(I), EPs_max_local)
310              ENDDO   ! end do (i=1,smax)
311     
312     ! for the case of all phases having same diameter
313     
314              IF (EPs_max_local == ONE) EPs_max_local = EPs_max_TMP(1)
315              CALC_EP_star = ONE - EPs_max_local
316     ! end YU_STANDISH section
317     ! ----------------------------------------------------------------<<<
318     
319     ! Part implemented by Dinesh for binary mixture, uncomment to use (Sof)
320     
321     !       if ((EP_s(IJK,1)+EP_s(IJK,2)) .NE. ZERO) THEN
322     !          xbar = EP_s(IJK,1)/(EP_s(IJK,1)+EP_s(IJK,2))
323     
324     !          if (xbar .LE. ep_s_max_ratio(1,2)) THEN
325     !             CALC_EP_star =MAX(0.36d0, (ONE-(((ep_s_max(1)-ep_s_max(2))+&
326     !              (ONE-d_p_ratio(1,2))*(ONE-ep_s_max(1))*ep_s_max(2))*(ep_s_max(1)+&
327     !              (ONE-ep_s_max(1)) *ep_s_max(2))*xbar/ep_s_max(1)+ep_s_max(2))))
328     !          else
329     !             CALC_EP_star =MAX(0.36d0, (ONE-((ONE -d_p_ratio(1,2))*(ep_s_max(1)&
330     !              +(ONE-ep_s_max(1))*ep_s_max(2))*(ONE -xbar) +ep_s_max(1))))
331     !          end if
332     !       else
333     !          CALC_EP_star = ONE - MIN(ep_s_max(1), ep_s_max(2)) !corrected by sof
334     !       end if
335     
336     ! Use the code (below) instead of the above commented code because the
337     ! phases were not rearranged and I didn't want to modify it (sof)
338     ! If you don't understand what's going on, contact me: sof@fluent.com
339     
340     ! In the case of binary mixture (Fedors-Landel empirical correlation)
341     ! ---------------------------------------------------------------->>>
342           ELSEIF(FEDORS_LANDEL) THEN
343     
344              IF(COMP_X_I(1) .LE. (EPs_max_TMP(1)/(EPs_max_TMP(1)+ &
345                                  (ONE - EPs_max_TMP(1))*EPs_max_TMP(2))) ) THEN
346     
347                 CALC_EP_star = (EPs_max_TMP(1) - EPs_max_TMP(2) + &
348                    (1 - sqrt(R_IJ(2,1))) * (ONE - EPs_max_TMP(1)) * &
349                    EPs_max_TMP(2) )*&
350                    (EPs_max_TMP(1) + (ONE - EPs_max_TMP(1)) * &
351                    EPs_max_TMP(2)) * COMP_X_I(1)/EPs_max_TMP(1) + &
352                    EPs_max_TMP(2)
353              ELSE
354                 CALC_EP_star = (ONE-sqrt(R_IJ(2,1))) * (EPs_max_TMP(1)+&
355                    (ONE-EPs_max_TMP(1)) * EPs_max_TMP(2)) * &
356                    (ONE - COMP_X_I(1)) + EPs_max_TMP(1)
357              ENDIF
358     ! this is gas volume fraction at packing
359              CALC_EP_star = ONE - CALC_EP_star
360            ENDIF ! for Yu_Standish and Fedors_Landel correlations
361     ! end FEDORS_LANDEL correlation
362     ! ----------------------------------------------------------------<<<
363     
364           RETURN
365           END FUNCTION CALC_ep_star
366