File: /nfs/home/0/users/jenkins/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           DOUBLE PRECISION, EXTERNAL :: BLEND_FUNCTION
70     !-----------------------------------------------
71     
72     !!$omp parallel do private(ijk)
73     !!   HPF$ independent
74     
75           DO IJK = ijkstart3, ijkend3
76              IF (FLUID_AT(IJK)) THEN
77     
78     
79                 IF (YU_STANDISH .OR. FEDORS_LANDEL) THEN
80     ! if Yu_Standish or Fedors_Landel correlations are used, then
81     ! ep_star_array is modified. this is the only time ep_star_array is
82     ! modified (see set_constprop).  (sof Nov-16-2005)
83                    EP_star_array(ijk) = calc_ep_star(ijk)
84     
85     ! now the values of ep_g_blend_start and ep_g_blend_end need to be
86     ! reassigned based on the new values of ep_star_array
87                    IF(BLENDING_STRESS.AND.TANH_BLEND) THEN
88                       ep_g_blend_start(ijk) = ep_star_array(ijk) * 0.99d0
89                       ep_g_blend_end(ijk)   = ep_star_array(ijk) * 1.01d0
90                    ELSEIF(BLENDING_STRESS.AND.SIGM_BLEND) THEN
91                       ep_g_blend_start(ijk) = EP_star_array(ijk) * 0.97d0
92                       ep_g_blend_end(ijk) = EP_star_array(ijk) * 1.01d0
93                    ELSE
94                       ep_g_blend_start(ijk) = ep_star_array(ijk)
95                       ep_g_blend_end(ijk)   = ep_star_array(ijk)
96                    ENDIF
97                 ENDIF
98     
99                 IF (EP_G(IJK) < EP_g_blend_end(ijk)) THEN
100                    P_STAR(IJK) = NEG_H(EP_G(IJK),EP_g_blend_end(ijk))
101                    IF(BLENDING_STRESS) THEN
102                       blend =  blend_function(IJK)
103                       P_STAR (IJK) = (1.0d0-blend) * P_STAR (IJK)
104                    ENDIF
105                 ELSE
106                    P_STAR(IJK) = ZERO
107                 ENDIF
108              ENDIF
109           ENDDO
110     
111           RETURN
112           END SUBROUTINE CALC_P_STAR
113     
114     
115     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
116     !                                                                      C
117     !  FUNCTION: CALC_ep_star                                              C
118     !  Purpose: calculate the local value of maximum packing               C
119     !                                                                      C
120     !  Author: D. Gera and M. Syamlal                     Date: 31-DEC-02  C
121     !  Reviewer:                                          Date:            C
122     !  Modified: S. Benyahia                              Date: 02-May-05  C
123     !                                                                      C
124     !  Literature/Document References:                                     C
125     !    A.B. Yu and N. Standish. Powder Tech, 52 (1987) 233-241           C
126     !    R.F. Fedors and R.F. Landel. Powder Tech, 23 (1979) 225-231       C
127     !                                                                      C
128     !  Variables referenced:                                               C
129     !  Variables modified:                                                 C
130     !  Local variables:                                                    C
131     !                                                                      C
132     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
133     
134           DOUBLE PRECISION FUNCTION CALC_ep_star(IJK)
135     
136     !-----------------------------------------------
137     ! Modules
138     !-----------------------------------------------
139           USE param
140           USE param1
141           USE fldvar
142           USE geometry
143           USE indices
144           USE physprop
145           USE constant
146           USE toleranc
147           USE compar
148           USE run
149           USE fun_avg
150           USE functions
151           IMPLICIT NONE
152     !-----------------------------------------------
153     ! Dummy arguments
154     !-----------------------------------------------
155     ! IJK index
156           INTEGER, INTENT(IN) :: IJK
157     !-----------------------------------------------
158     ! Local variables
159     !-----------------------------------------------
160     ! Indices
161           INTEGER :: I, J
162     
163     !      DOUBLE PRECISION :: xbar
164     
165     ! start sof modifications (02-May-05)
166     ! maximum packing for the mixture
167            DOUBLE PRECISION :: P_IT(MMAX)
168     ! true maximum packing for the mixture
169            DOUBLE PRECISION :: EPs_max_local
170     ! maximum packing fraction for a binary mixture
171            DOUBLE PRECISION :: P_IJ(MMAX, MMAX)
172     ! particle diameter ratio
173            DOUBLE PRECISION :: R_IJ(MMAX, MMAX)
174     ! fractional solids volume corresponding to P_IJ
175            DOUBLE PRECISION :: X_IJ(MMAX, MMAX)
176     ! fractional solids volume in a mixture
177     ! this is Xj in eq. 22 of Yu-Standish
178            DOUBLE PRECISION :: COMP_X_I(MMAX), SUM_LOCAL
179     ! local aliases for particle diameter, solids volume fraction and the
180     ! maximum solids volume fraction which are used to rearrange solids
181     ! phases from coarsest to finest
182            DOUBLE PRECISION :: DP_TMP(MMAX), EPs_TMP(MMAX), &
183                                EPs_max_TMP(MMAX), old_value
184     !-----------------------------------------------
185     
186     
187           IF (CALL_DQMOM) THEN
188     ! sort particles to start from coarsest to finest particles
189     ! assigning values to local aliases
190              DO I = 1, SMAX
191                 DP_TMP(I) = D_P(IJK,I)
192                 EPs_TMP(I) = EP_s(IJK,I)
193                 EPs_max_TMP(I) = ep_s_max(I)
194              ENDDO
195     ! sorting particles from coarse to fine
196              DO I = 1, SMAX
197                 DO J = I , SMAX
198     ! check if phase J is larger than phase I
199                    IF(DP_TMP(I) < DP_TMP(J)) THEN
200     ! temporarily store phase i diameter
201                       old_value = DP_TMP(I)
202     ! overwrite phase i diameter with smaller phase j diameter
203                       DP_TMP(I) = DP_TMP(J)
204     ! overwrite phase j diameter with stoired phase i diameter
205                       DP_TMP(J) = old_value
206     
207                       old_value = EPs_TMP(I)
208                       EPs_TMP(I) = EPs_TMP(J)
209                       EPs_TMP(J) = old_value
210     
211                       old_value = EPs_max_TMP(I)
212                       EPs_max_TMP(I) = EPs_max_TMP(J)
213                       EPs_max_TMP(J) = old_value
214                    ENDIF
215                 ENDDO
216              ENDDO
217     
218           ELSE  ! not dqmom
219     
220     ! assigning values to local aliases
221              DO I = 1, SMAX
222                 DP_TMP(I) = D_P(IJK,M_MAX(I))
223                 EPs_TMP(I) = EP_s(IJK,M_MAX(I))
224                 EPs_max_TMP(I) = ep_s_max(M_MAX(I))
225              ENDDO
226           ENDIF   ! end if/else (call_dqmom)
227     
228     
229     ! this is the way the algorithm was written by Yu and Standish (sof).
230     ! compute equations 25 in Yu-Standish
231     ! (this is also needed by Fedors_Landel)
232           DO I = 1, SMAX
233              SUM_LOCAL = ZERO
234              DO J = 1, SMAX
235                 IF(I .GE. J) THEN
236                    R_IJ(I,J) = DP_TMP(I)/DP_TMP(J)
237                 ELSE
238                    R_IJ(I,J) = DP_TMP(J)/DP_TMP(I)
239                 ENDIF
240                 SUM_LOCAL = SUM_LOCAL + EPs_TMP(J)
241              ENDDO   ! end do (j=1,smax)
242     
243              IF(SUM_LOCAL > DIL_EP_s) THEN
244     ! fractional solids volume see eq. 20
245                 COMP_X_I(I) = EPs_TMP(I)/SUM_LOCAL
246              ELSE
247     ! return first phase ep_s_max in case very dilute
248                 CALC_EP_star = ONE - EPs_max_TMP(1)
249                 RETURN
250              ENDIF
251           ENDDO   ! end do (i=1,smax)
252     
253     ! Begin YU_STANDISH section
254     ! ---------------------------------------------------------------->>>
255           IF(YU_STANDISH) THEN
256     ! compute equation 23-24 in Yu-Standish
257              DO I = 1, SMAX
258                 DO J = 1, SMAX
259                    IF(R_IJ(I,J) .LE. 0.741d0) THEN
260                       IF(J .LT. I) THEN
261                          X_IJ(I,J) = (ONE - R_IJ(I,J)*R_IJ(I,J))/&
262                                      (2.0d0 -  EPs_max_TMP(I))
263                       ELSE
264                          X_IJ(I,J) = ONE - (ONE - R_IJ(I,J)*R_IJ(I,J))/&
265                                      (2.0d0 -  EPs_max_TMP(I))
266                      ENDIF
267                      P_IJ(I, J) = EPs_max_TMP(I) + EPs_max_TMP(I)*&
268                         (ONE-EPs_max_TMP(I)) * (ONE - 2.35d0*R_IJ(I,J) + &
269                         1.35d0*R_IJ(I,J)*R_IJ(I,J) )
270                    ELSE
271                       P_IJ(I, J) = EPs_max_TMP(I)
272                    ENDIF
273                 ENDDO   ! end do (j=1,smax)
274              ENDDO   ! end do (i=1,smax)
275     
276     ! Compute equation 22
277              EPs_max_local = ONE
278              DO I = 1, SMAX
279                 SUM_LOCAL = ZERO
280     
281                 IF(I .GE. 2) THEN
282                    DO J = 1, (I-1)
283                       IF(P_IJ(I,J) == EPs_max_TMP(I)) THEN
284                          SUM_LOCAL = SUM_LOCAL
285                       ELSE
286                          SUM_LOCAL = SUM_LOCAL + (ONE - EPs_max_TMP(I)/&
287                             P_IJ(I,J))*COMP_X_I(J)/X_IJ(I,J)
288                       ENDIF
289                    ENDDO
290                 ENDIF
291     
292                 IF((I+1) .LE. SMAX) THEN
293                    DO J = (I+1), SMAX
294                       IF( P_IJ(I, J) == EPs_max_TMP(I) ) THEN
295                          SUM_LOCAL = SUM_LOCAL
296                       ELSE
297                          SUM_LOCAL = SUM_LOCAL + (ONE - EPs_max_TMP(I)/&
298                             P_IJ(I, J))*COMP_X_I(J)/X_IJ(I, J)
299                       ENDIF
300                    ENDDO
301                 ENDIF
302     
303                 IF (SUM_LOCAL .NE. ZERO) THEN
304                    P_IT(I) = EPs_max_TMP(I)/(ONE - SUM_LOCAL)
305                 ELSE
306     ! do nothing if particles have same diameter
307                    P_IT(I) = ONE
308                 ENDIF
309     
310                 EPs_max_local = MIN(P_IT(I), EPs_max_local)
311              ENDDO   ! end do (i=1,smax)
312     
313     ! for the case of all phases having same diameter
314     
315              IF (EPs_max_local == ONE) EPs_max_local = EPs_max_TMP(1)
316              CALC_EP_star = ONE - EPs_max_local
317     ! end YU_STANDISH section
318     ! ----------------------------------------------------------------<<<
319     
320     ! Part implemented by Dinesh for binary mixture, uncomment to use (Sof)
321     
322     !       if ((EP_s(IJK,1)+EP_s(IJK,2)) .NE. ZERO) THEN
323     !          xbar = EP_s(IJK,1)/(EP_s(IJK,1)+EP_s(IJK,2))
324     
325     !          if (xbar .LE. ep_s_max_ratio(1,2)) THEN
326     !             CALC_EP_star =MAX(0.36d0, (ONE-(((ep_s_max(1)-ep_s_max(2))+&
327     !              (ONE-d_p_ratio(1,2))*(ONE-ep_s_max(1))*ep_s_max(2))*(ep_s_max(1)+&
328     !              (ONE-ep_s_max(1)) *ep_s_max(2))*xbar/ep_s_max(1)+ep_s_max(2))))
329     !          else
330     !             CALC_EP_star =MAX(0.36d0, (ONE-((ONE -d_p_ratio(1,2))*(ep_s_max(1)&
331     !              +(ONE-ep_s_max(1))*ep_s_max(2))*(ONE -xbar) +ep_s_max(1))))
332     !          end if
333     !       else
334     !          CALC_EP_star = ONE - MIN(ep_s_max(1), ep_s_max(2)) !corrected by sof
335     !       end if
336     
337     ! Use the code (below) instead of the above commented code because the
338     ! phases were not rearranged and I didn't want to modify it (sof)
339     ! If you don't understand what's going on, contact me: sof@fluent.com
340     
341     ! In the case of binary mixture (Fedors-Landel empirical correlation)
342     ! ---------------------------------------------------------------->>>
343           ELSEIF(FEDORS_LANDEL) THEN
344     
345              IF(COMP_X_I(1) .LE. (EPs_max_TMP(1)/(EPs_max_TMP(1)+ &
346                                  (ONE - EPs_max_TMP(1))*EPs_max_TMP(2))) ) THEN
347     
348                 CALC_EP_star = (EPs_max_TMP(1) - EPs_max_TMP(2) + &
349                    (1 - sqrt(R_IJ(2,1))) * (ONE - EPs_max_TMP(1)) * &
350                    EPs_max_TMP(2) )*&
351                    (EPs_max_TMP(1) + (ONE - EPs_max_TMP(1)) * &
352                    EPs_max_TMP(2)) * COMP_X_I(1)/EPs_max_TMP(1) + &
353                    EPs_max_TMP(2)
354              ELSE
355                 CALC_EP_star = (ONE-sqrt(R_IJ(2,1))) * (EPs_max_TMP(1)+&
356                    (ONE-EPs_max_TMP(1)) * EPs_max_TMP(2)) * &
357                    (ONE - COMP_X_I(1)) + EPs_max_TMP(1)
358              ENDIF
359     ! this is gas volume fraction at packing
360              CALC_EP_star = ONE - CALC_EP_star
361            ENDIF ! for Yu_Standish and Fedors_Landel correlations
362     ! end FEDORS_LANDEL correlation
363     ! ----------------------------------------------------------------<<<
364     
365           RETURN
366           END FUNCTION CALC_ep_star
367     
368     
369     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
370     !                                                                      C
371     !  FUNCTION: blend_function                                            C
372     !  Purpose: To calculate blending function                             C
373     !                                                                      C
374     !  Author: S. Pannala                                 Date: 28-FEB-06  C
375     !  Reviewer:                                          Date:            C
376     !  Modified:                                          Date: 24-OCT-06  C
377     !                                                                      C
378     !                                                                      C
379     !  Literature/Document References:                                     C
380     !  Variables referenced:                                               C
381     !  Variables modified:                                                 C
382     !                                                                      C
383     !  Local variables:                                                    C
384     !                                                                      C
385     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
386     
387           DOUBLE PRECISION FUNCTION blend_function(IJK)
388     
389     !-----------------------------------------------
390     ! Modules
391     !-----------------------------------------------
392           USE param
393           USE param1
394           USE fldvar
395           USE geometry
396           USE indices
397           USE physprop
398           USE constant
399           USE toleranc
400           USE compar
401           USE run
402           USE visc_s
403           USE fun_avg
404           USE functions
405           IMPLICIT NONE
406     !-----------------------------------------------
407     ! Dummy arguments
408     !-----------------------------------------------
409     ! IJK index
410           INTEGER, INTENT(IN) :: IJK
411     !-----------------------------------------------
412     ! Local variables
413     !-----------------------------------------------
414     ! Logical to see whether this is the first entry to this routine
415           LOGICAL,SAVE:: FIRST_PASS = .TRUE.
416     ! Blend Factor
417           Double Precision:: blend, blend_right
418     ! Scale Factor
419           Double Precision, Save:: scale
420     ! Midpoint
421           Double Precision, Save:: ep_mid_point
422     !-----------------------------------------------
423     
424     ! Tan hyperbolic blending of stresses
425           IF(TANH_BLEND) THEN
426              IF(EP_g(IJK) .LT. ep_g_blend_end(ijk).AND. EP_g(IJK) .GT. ep_g_blend_start(ijk)) THEN
427                 ep_mid_point = (ep_g_blend_end(IJK)+ep_g_blend_start(IJK))/2.0d0
428                 blend = tanh(2.0d0*pi*(ep_g(IJK)-ep_mid_point)/ &
429                 (ep_g_blend_end(IJK)-ep_g_blend_start(IJK)))
430                 blend = (blend+1.0d0)/2.0d0
431              ELSEIF(EP_g(IJK) .GE. ep_g_blend_end(ijk)) THEN
432                 blend = 1.0d0
433              ELSEIF(EP_g(IJK) .LE. ep_g_blend_start(ijk)) THEN
434                 blend = 0.0d0
435              ENDIF
436     
437     ! Truncated and Scaled Sigmoidal blending of stresses
438           ELSEIF(SIGM_BLEND) THEN
439              IF(FIRST_PASS) THEN
440                 blend_right =  1.0d0/(1+0.01d0**((ep_g_blend_end(IJK)-ep_star_array(IJK))&
441                 /(ep_g_blend_end(IJK)-ep_g_blend_start(IJK))))
442                 blend_right = (blend_right+1.0d0)/2.0d0
443                 scale = 1.0d0/blend_right
444                 write(*,*) 'Blending value at end and scaling factor', blend_right, scale
445                 FIRST_PASS = .FALSE.
446              ENDIF
447              IF(EP_g(IJK) .LT. ep_g_blend_end(ijk)) THEN
448                 blend =  scale/(1+0.01d0**((ep_g(IJK)-ep_star_array(IJK))&
449                 /(ep_g_blend_end(IJK)-ep_g_blend_start(IJK))))
450              ELSEIF(EP_g(IJK) .GE. ep_g_blend_end(ijk)) THEN
451                 blend = 1.0d0
452              ENDIF
453     
454           ENDIF
455     
456           blend_function = blend
457     
458           RETURN
459           END FUNCTION blend_function
460     
461