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

1     ! TO DO:
2     ! 1. Check the formulation based on MCP.
3     ! 2. The pressure correction should be based on sum of close-packed
4     !    solids?
5     
6     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
7     !                                                                      C
8     !  Subroutine: CALC_VOL_FR                                             C
9     !  Purpose: Calculate volume fractions of phases used in pressure      C
10     !           corrections.                                               C
11     !                                                                      C
12     !  Notes: see mark_phase_4_cor for more details                        C
13     !                                                                      C
14     !  Author: M. Syamlal                                 Date: 5-JUL-96   C
15     !                                                                      C
16     !                                                                      C
17     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
18           SUBROUTINE CALC_VOL_FR(P_STAR, RO_G, ROP_G, EP_G, ROP_S, IER)
19     
20     
21     ! Modules
22     !---------------------------------------------------------------------//
23           USE param
24           USE param1
25           USE parallel
26           USE run
27           USE geometry
28           USE indices
29           USE physprop
30           USE visc_s
31           USE constant
32           USE pgcor
33           USE pscor
34           USE compar
35           USE sendrecv
36           USE discretelement
37           USE mpi_utility
38           use fldvar, only: RO_S, EP_S
39           USE solids_pressure
40           USE functions
41     
42           IMPLICIT NONE
43     
44     ! Dummy Arguments
45     !---------------------------------------------------------------------//
46     ! Solids pressure
47           DOUBLE PRECISION, INTENT(IN) :: P_star(DIMENSION_3)
48     ! Gas density
49           DOUBLE PRECISION, INTENT(INOUT) :: RO_g(DIMENSION_3)
50     ! Gas bulk density
51           DOUBLE PRECISION, INTENT(INOUT) :: ROP_g(DIMENSION_3)
52     ! Gas volume fraction
53           DOUBLE PRECISION, INTENT(INOUT) :: EP_g(DIMENSION_3)
54     ! solids bulk densities
55           DOUBLE PRECISION, INTENT(INOUT) :: ROP_s(DIMENSION_3, DIMENSION_M)
56     ! error index
57           INTEGER, INTENT(INOUT) :: IER
58     
59     ! Local variables
60     !---------------------------------------------------------------------//
61     ! volume of particle type M for GHD theory
62           DOUBLE PRECISION :: VOL_M
63     ! volume fraction of close-packed region
64           DOUBLE PRECISION :: EPcp
65     ! volume fraction of solids phase
66           DOUBLE PRECISION :: EPS
67     ! sum of volume fractions
68           DOUBLE PRECISION :: SUMVF
69     ! Whichever phase is given by MF becomes the phase whose volume fraction
70     ! is corrected based on all other phases present.  Generally MF gets
71     ! defined as 0 so that the gas phase void fraction becomes corrected
72     ! based on the summation of the volume fractions of all other phases
73           INTEGER :: MF
74     ! Whichever phase is given by MCPl becomes the phase whose volume
75     ! fraction is corrected based on value of maximum close packing and all
76     ! other solids phase that can close pack.  This is basically a local
77     ! variable for MCP but only in cells where close packed conditions
78     ! exist.
79           INTEGER :: MCPl
80     ! Index of solids phase
81           INTEGER :: M
82     ! Indices
83           INTEGER :: IJK
84     
85     ! Arrays for storing errors:
86     ! 110 - Negative volume fraciton
87     ! 11x - Unclassified
88           INTEGER :: Err_l(0:numPEs-1)  ! local
89           INTEGER :: Err_g(0:numPEs-1)  ! global
90     
91           LOGICAL, PARAMETER :: REPORT_NEG_VOLFRAC = .TRUE.
92     
93     ! Flag to write log header
94           LOGICAL :: wHeader
95     !---------------------------------------------------------------------//
96     
97     
98     ! Initialize:
99           wHeader = .TRUE.
100     ! Initialize error flag.
101           Err_l = 0
102     
103     !!$omp  parallel do private(MCPl, EPCP, SUMVF, MF, M) &
104     !!$omp&  schedule(static)
105     
106           DO IJK = ijkstart3, ijkend3
107              IF (FLUID_AT(IJK)) THEN
108     
109     ! calculate the volume fraction of the solids phase based on the volume
110     ! fractions of all other solids phases and on the value of maximum
111     ! packing when that solids phase continuity was not solved for the
112     ! indicated cell.
113     !----------------------------------------------------------------->>>
114                 IF (PHASE_4_P_S(IJK) /= UNDEFINED_I) THEN
115     ! for the current cell check if a solids phase continuity was skipped.
116     ! this value will either be undefined (no cell was skipped in any
117     ! solids continuity equations) or be a solids phase index (indicated
118     ! solids continuity was skipped) for a solids phase that does close pack
119     ! (i.e., close_packed=T) and the cell exhibits close packing conditions.
120     ! note: this branch is never entered given the existing version of
121     ! mark_phase_4_cor.
122                    MCPl   = PHASE_4_P_s(IJK)
123     
124     ! finding the value of maximum close packing based on the expression for
125     ! plastic pressure
126                    EPCP = 1. - INV_H(P_STAR(IJK),EP_g_blend_end(ijk))
127     ! summing the solids volume fraction of all continuum solids phases
128     ! that are marked as close_packed except for the solids phase which
129     ! is also marked by phase_4_p_s (skipping its own continuity)
130                    SUMVF = ZERO
131                    DO M = 1, MMAX
132                       IF (CLOSE_PACKED(M) .AND. M/=MCPl) SUMVF = SUMVF + EP_S(IJK,M)
133                    ENDDO
134     ! sum of solids volume fraction from the DEM
135                    IF (DES_CONTINUUM_HYBRID) THEN
136                       DO M = 1,DES_MMAX
137                          EPS = DES_ROP_S(IJK,M)/DES_RO_S(M)
138                          SUMVF = SUMVF + EPS
139                       ENDDO
140                    ENDIF
141                    ROP_S(IJK,MCPl) = (EPCP - SUMVF)*RO_S(IJK,MCPl)
142                 ENDIF
143     !-----------------------------------------------------------------<<<
144     
145     
146     ! calculate the volume fraction of the 'solids' phase based on the
147     ! volume fractions of all other phases (including gas) if the gas
148     ! continuity was solved rather than that phases own continuity.
149     ! if the gas continuity was not solved then calculate the void
150     ! fraction of the gas phase based on all other solids phases.
151     !----------------------------------------------------------------->>>
152     ! for the current cell check if the gas phase continuity was solved for
153     ! the gas phase while the solids phase continuity for the indicated
154     ! solids phase was skipped. this value will either be 0 (gas phase
155     ! continuity was not solved) or be a solids phase index (gas continuity
156     ! was solved while the indicated solids phase continuity was skipped)
157     ! for a solids phase that does not close pack (i.e., close_packed=F) and
158     ! is in greater concentration than the gas phase.
159     ! Note: MF will always be set to 0 here given the existing version of
160     ! mark_phase_4_cor
161                 MF = PHASE_4_P_G(IJK)
162     
163                 SUMVF = ZERO
164                 IF (0 /= MF) THEN
165     ! if gas continuity was solved rather than the solids phase, then
166     ! include the gas phase void fraction in the summation here.
167                    EP_G(IJK) = ROP_G(IJK)/RO_G(IJK)
168                    SUMVF = SUMVF + EP_G(IJK)
169                 ENDIF
170     
171     ! modified for GHD theory
172                 IF(KT_TYPE_ENUM == GHD_2007) THEN
173                   ROP_S(IJK,MMAX) = ZERO  ! mixture density
174                   DO M = 1, SMAX
175     ! volume of particle M based on fixed diamter Dp0
176                      VOL_M = PI*D_P0(M)**3/6d0
177                      IF (M /= MF) THEN
178                        SUMVF = SUMVF + EP_S(IJK,M)
179                        ROP_S(IJK,MMAX) = ROP_S(IJK,MMAX) + RO_S(IJK,M)*EP_S(IJK,M)
180                      ENDIF
181                   ENDDO
182                 ELSE
183     ! summing the solids volume fraction of all continuum solids phases
184     ! except for the continuum solids phase which was marked by phase_4_p_g
185     ! (skipping its continuity while solving gas continuity)
186                   DO M = 1, MMAX
187                      IF (M /= MF) SUMVF = SUMVF + EP_S(IJK,M)
188                   ENDDO
189                 ENDIF   ! end if/else trim(kt_type)=='ghd'
190     
191     ! sum of solids volume fraction from the DEM
192                 IF (DES_CONTINUUM_HYBRID) THEN
193                    DO M = 1,DES_MMAX
194                       EPS = DES_ROP_S(IJK,M)/DES_RO_S(M)
195                       SUMVF = SUMVF + EPS
196                    ENDDO
197                 ENDIF
198     
199                 IF (0 == MF) THEN
200     ! if no gas phase continuity was solved in the current cell then correct
201     ! the void fraction of the gas phase based on the total solids volume
202     ! fraction of all solids phases
203                    EP_G(IJK) = ONE - SUMVF
204     
205     ! Set error flag for negative volume fraction.
206                    IF (EP_G(IJK) < ZERO) THEN
207                       Err_l(myPE) = 110
208                       IF(REPORT_NEG_VOLFRAC) CALL EPgErr_LOG(IJK, wHeader)
209                    ENDIF
210     
211                    ROP_G(IJK) = EP_G(IJK)*RO_G(IJK)
212                 ELSE
213     ! else correct the volume fraction of the solids phase that was marked
214                    ROP_S(IJK,MF) = (ONE - SUMVF)*RO_S(IJK,MF)
215                 ENDIF
216     !-----------------------------------------------------------------<<<
217     
218              ENDIF    ! end if (fluid_at(ijk))
219           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
220     
221           CALL send_recv(EP_G, 2)
222           CALL send_recv(ROP_G, 2)
223           CALL send_recv(ROP_S, 2)
224     
225           CALL global_all_sum(Err_l, Err_g)
226           IER = maxval(Err_g)
227     
228     
229           RETURN
230           END SUBROUTINE CALC_VOL_FR
231     
232     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
233     !                                                                      C
234     !  Subroutine: SET_EP_FACTORS                                          C
235     !  Purpose: Set multiplication factors to ep_g to control the form     C
236     !  of the governing equations                                          C
237     !                                                                      C
238     !  References:                                                         C
239     !  Anderson and Jackson, A fluid mechanical description of fluidized   C
240     !     beds, Ind. Eng. Chem. Fundam., 6(4) 1967, 527-539.               C
241     !  Ishii, Thermo-fluid dynamic theory of two-phase flow, des Etudes    C
242     !     et Recherches D'Electricite de France, Eyrolles, Paris, France   C
243     !     1975.                                                            C
244     !                                                                      C
245     !                                                                      C
246     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
247           SUBROUTINE SET_EP_FACTORS
248     
249     ! Modules
250     !---------------------------------------------------------------------//
251           use compar, only: ijkstart3, ijkend3
252           use fldvar, only: ep_g, ep_s
253           use fldvar, only: epg_ifac, eps_ifac, epg_jfac
254           use param1, only: one, undefined
255           use physprop, only: mmax, mu_s0
256           use run, only: jackson, ishii
257           IMPLICIT NONE
258     ! Local variables
259     !---------------------------------------------------------------------//
260     ! indices
261           INTEGER :: ijk, m
262     !---------------------------------------------------------------------//
263           epg_jfac(:) = one
264           epg_ifac(:) = one
265           eps_ifac(:,:) = one
266     
267           if (jackson) then
268              do ijk = ijkstart3, ijkend3
269                 epg_jfac(ijk) = ep_g(ijk)
270              enddo
271           endif
272           if (ishii) then
273              do ijk = ijkstart3, ijkend3
274                 epg_ifac(ijk) = ep_g(ijk)
275                 do m = 1, mmax
276     ! This seems more appropriate for non-granular systems
277     ! So currently only incorporate this for 'constant' viscosity cases
278                    IF (mu_s0(m) /= undefined) THEN
279                       eps_ifac(ijk,m) = ep_s(ijk,m)
280                    ENDIF
281                 enddo
282              enddo
283           endif
284     
285           RETURN
286           END SUBROUTINE SET_EP_FACTORS
287     
288     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
289     !                                                                      !
290     !  Subroutine: EPgErr_LOG                                              !
291     !  Purpose: Record information about the location and conditions that  !
292     !           resulted in a negative gas phase volume fraction.          !
293     !                                                                      !
294     !  Author: J. Musser                                  Date: 16-APR-15  !
295     !  Reviewer:                                          Date:            !
296     !                                                                      !
297     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
298           SUBROUTINE EPgErr_LOG(IJK, tHeader)
299     
300     ! Simulation time
301           use run, only: TIME
302     ! Gas phase pressure.
303           use fldvar, only: EP_g, EP_s, ROP_s, RO_s, U_s, V_s, W_s
304           use cutcell
305     
306           use physprop, only: SMAX
307           use compar
308           use indices
309           use functions
310           use error_manager
311     
312     
313           IMPLICIT NONE
314     
315           INTEGER, intent(in) :: IJK
316           LOGICAL, intent(inout) :: tHeader
317     
318           LOGICAL :: lExists
319           CHARACTER(LEN=255) :: lFile
320           INTEGER, parameter :: lUnit = 4868
321           LOGICAL, save :: fHeader = .TRUE.
322     
323           INTEGER :: M
324     
325           lFile = '';
326           if(numPEs > 1) then
327              write(lFile,"('EPgErr_',I4.4,'.log')") myPE
328           else
329              write(lFile,"('EPgErr.log')")
330           endif
331           inquire(file=trim(lFile),exist=lExists)
332           if(lExists) then
333              open(lUnit,file=trim(adjustl(lFile)),                         &
334                 status='old', position='append')
335           else
336              open(lUnit,file=trim(adjustl(lFile)), status='new')
337           endif
338     
339           if(fHeader) then
340              write(lUnit,1000)
341              fHeader = .FALSE.
342           endif
343     
344           if(tHeader) then
345              write(lUnit,"(/2x,'Simulation time: ',g12.5)") TIME
346              tHeader = .FALSE.
347           endif
348     
349           write(lUnit,1001) IJK, I_OF(IJK), J_OF(IJK), K_OF(IJK)
350           write(lUnit,"(6x,A,1X,g12.5)",ADVANCE='NO') 'EP_g:', EP_g(IJK)
351           DO M=1,SMAX
352              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
353                 trim(iVar('EP_s',M)), EP_s(IJK,M)
354           ENDDO
355           write(lUnit,"(' ')")
356     
357           write(lUnit,"(24x)", ADVANCE='NO')
358           DO M=1,SMAX
359              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
360                 trim(iVar('RO_s',M)), RO_s(IJK,M)
361           ENDDO
362           write(lUnit,"(' ')")
363     
364           write(lUnit,"(24x)", ADVANCE='NO')
365           DO M=1,SMAX
366              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
367                 trim(iVar('ROPs',M)), ROP_s(IJK,M)
368           ENDDO
369           write(lUnit,"(24x,' ')")
370     
371     
372           write(lUnit,"(24x)", ADVANCE='NO')
373           DO M=1,SMAX
374              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
375                 trim(iVar('U_se',M)), U_s(IJK,M)
376           ENDDO
377           write(lUnit,"(24x,' ')")
378     
379           write(lUnit,"(24x)", ADVANCE='NO')
380           DO M=1,SMAX
381              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
382                 trim(iVar('V_sn',M)), V_s(IJK,M)
383           ENDDO
384           write(lUnit,"(24x,' ')")
385     
386           write(lUnit,"(24x)", ADVANCE='NO')
387           DO M=1,SMAX
388              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
389                 trim(iVar('W_st',M)), W_s(IJK,M)
390           ENDDO
391           write(lUnit,"(24x,' ')")
392     
393           write(lUnit,"(24x)", ADVANCE='NO')
394           DO M=1,SMAX
395              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
396                 trim(iVar('U_sw',M)), U_s(WEST_OF(IJK),M)
397           ENDDO
398           write(lUnit,"(24x,' ')")
399     
400           write(lUnit,"(24x)", ADVANCE='NO')
401           DO M=1,SMAX
402              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
403                 trim(iVar('V_ss',M)), V_s(SOUTH_OF(IJK),M)
404           ENDDO
405           write(lUnit,"(24x,' ')")
406     
407           write(lUnit,"(24x)", ADVANCE='NO')
408           DO M=1,SMAX
409              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
410                 trim(iVar('W_sb',M)), W_s(BOTTOM_OF(IJK),M)
411           ENDDO
412           write(lUnit,"(24x,' ')")
413     
414     
415           if(CARTESIAN_GRID) then
416              write(lUnit,"(6x,A,1X,L1)",ADVANCE='NO') 'Cut Cell:', CUT_CELL_AT(IJK)
417              write(lUnit,"(6x,A,1X,L1)") 'Small Cell:', SMALL_CELL_AT(IJK)
418              write(lUnit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
419                 xg_e(I_OF(IJK)), yg_n(J_of(ijk)), zg_t(k_of(ijk))
420           endif
421     
422           close(lUnit)
423     
424           RETURN
425     
426      1000 FORMAT('One or more cells have reported a negative gas volume ', &
427              'fraction (EP_g).',/)
428     
429      1001 FORMAT(/4X,'IJK: ',I8,7X,'I: ',I4,'  J: ',I4,'  K: ',I4)
430     
431           END SUBROUTINE EPgErr_LOG
432