File: N:\mfix\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 compar
24           USE constant
25           USE discretelement
26           USE functions
27           USE geometry
28           USE indices
29           USE mpi_utility
30           USE parallel
31           USE param
32           USE param1
33           USE pgcor
34           USE physprop
35           USE pscor
36           USE run, only: ghd_2007, kt_type_enum
37           USE sendrecv
38           USE solids_pressure
39           USE visc_s
40           use fldvar, only: RO_S, EP_S
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+DES_MMAX
132                       IF (CLOSE_PACKED(M) .AND. M/=MCPl) SUMVF = SUMVF + EP_S(IJK,M)
133                    ENDDO
134                    ROP_S(IJK,MCPl) = (EPCP - SUMVF)*RO_S(IJK,MCPl)
135                 ENDIF
136     !-----------------------------------------------------------------<<<
137     
138     
139     ! calculate the volume fraction of the 'solids' phase based on the
140     ! volume fractions of all other phases (including gas) if the gas
141     ! continuity was solved rather than that phases own continuity.
142     ! if the gas continuity was not solved then calculate the void
143     ! fraction of the gas phase based on all other solids phases.
144     !----------------------------------------------------------------->>>
145     ! for the current cell check if the gas phase continuity was solved for
146     ! the gas phase while the solids phase continuity for the indicated
147     ! solids phase was skipped. this value will either be 0 (gas phase
148     ! continuity was not solved) or be a solids phase index (gas continuity
149     ! was solved while the indicated solids phase continuity was skipped)
150     ! for a solids phase that does not close pack (i.e., close_packed=F) and
151     ! is in greater concentration than the gas phase.
152     ! Note: MF will always be set to 0 here given the existing version of
153     ! mark_phase_4_cor
154                 MF = PHASE_4_P_G(IJK)
155     
156                 SUMVF = ZERO
157                 IF (0 /= MF) THEN
158     ! if gas continuity was solved rather than the solids phase, then
159     ! include the gas phase void fraction in the summation here.
160                    EP_G(IJK) = ROP_G(IJK)/RO_G(IJK)
161                    SUMVF = SUMVF + EP_G(IJK)
162                 ENDIF
163     
164     ! evaluate mixture bulk density rop_s(mmax) for GHD theory
165                 IF(KT_TYPE_ENUM == GHD_2007) THEN
166                   ROP_S(IJK,MMAX) = ZERO  ! mixture density
167                   DO M = 1, SMAX
168                      IF (M /= MF) THEN
169                        ROP_S(IJK,MMAX) = ROP_S(IJK,MMAX) + RO_S(IJK,M)*EP_S(IJK,M)
170                      ENDIF
171                   ENDDO
172                 ENDIF
173     
174     ! summing the solids volume fraction of all continuum solids phases
175     ! except for the continuum solids phase which was marked by phase_4_p_g
176     ! (skipping its continuity while solving gas continuity)
177                 DO M = 1, DES_MMAX+MMAX
178                    IF(KT_TYPE_ENUM == GHD_2007 .AND. M == MMAX) CYCLE
179                    IF (M /= MF) SUMVF = SUMVF + EP_S(IJK,M)
180                 ENDDO
181     
182                 IF (0 == MF) THEN
183     ! if no gas phase continuity was solved in the current cell then correct
184     ! the void fraction of the gas phase based on the total solids volume
185     ! fraction of all solids phases
186                    EP_G(IJK) = ONE - SUMVF
187     
188     ! Set error flag for negative volume fraction.
189                    IF (EP_G(IJK) < ZERO) THEN
190                       Err_l(myPE) = 110
191                       IF(REPORT_NEG_VOLFRAC) CALL EPgErr_LOG(IJK, wHeader)
192                    ENDIF
193     
194                    ROP_G(IJK) = EP_G(IJK)*RO_G(IJK)
195                 ELSE
196     ! else correct the volume fraction of the solids phase that was marked
197                    ROP_S(IJK,MF) = (ONE - SUMVF)*RO_S(IJK,MF)
198                 ENDIF
199     !-----------------------------------------------------------------<<<
200     
201              ENDIF    ! end if (fluid_at(ijk))
202           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
203     
204           CALL send_recv(EP_G, 2)
205           CALL send_recv(ROP_G, 2)
206           CALL send_recv(ROP_S, 2)
207     
208           CALL global_all_sum(Err_l, Err_g)
209           IER = maxval(Err_g)
210     
211     
212           RETURN
213           END SUBROUTINE CALC_VOL_FR
214     
215     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
216     !                                                                      C
217     !  Subroutine: SET_EP_FACTORS                                          C
218     !  Purpose: Set multiplication factors to ep_g to control the form     C
219     !  of the governing equations                                          C
220     !                                                                      C
221     !  References:                                                         C
222     !  Anderson and Jackson, A fluid mechanical description of fluidized   C
223     !     beds, Ind. Eng. Chem. Fundam., 6(4) 1967, 527-539.               C
224     !  Ishii, Thermo-fluid dynamic theory of two-phase flow, des Etudes    C
225     !     et Recherches D'Electricite de France, Eyrolles, Paris, France   C
226     !     1975.                                                            C
227     !                                                                      C
228     !                                                                      C
229     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
230           SUBROUTINE SET_EP_FACTORS
231     
232     ! Modules
233     !---------------------------------------------------------------------//
234           use compar, only: ijkstart3, ijkend3
235           use fldvar, only: ep_g, ep_s
236           use fldvar, only: epg_ifac, eps_ifac, epg_jfac
237           use functions, only: wall_at
238           use param1, only: one, undefined
239           use physprop, only: mmax, mu_s0
240           use run, only: jackson, ishii
241           IMPLICIT NONE
242     ! Local variables
243     !---------------------------------------------------------------------//
244     ! indices
245           INTEGER :: ijk, m
246     !---------------------------------------------------------------------//
247           epg_jfac(:) = one
248           epg_ifac(:) = one
249           eps_ifac(:,:) = one
250     
251           if (jackson) then
252              do ijk = ijkstart3, ijkend3
253                 if (wall_at(ijk)) cycle
254                 epg_jfac(ijk) = ep_g(ijk)
255              enddo
256           endif
257           if (ishii) then
258              do ijk = ijkstart3, ijkend3
259                 if (wall_at(ijk)) cycle
260                 epg_ifac(ijk) = ep_g(ijk)
261                 do m = 1, mmax
262     ! This seems more appropriate for non-granular systems
263     ! So currently only incorporate this for 'constant' viscosity cases
264                    IF (mu_s0(m) /= undefined) THEN
265                       eps_ifac(ijk,m) = ep_s(ijk,m)
266                    ENDIF
267                 enddo
268              enddo
269           endif
270     
271           RETURN
272           END SUBROUTINE SET_EP_FACTORS
273     
274     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
275     !                                                                      !
276     !  Subroutine: EPgErr_LOG                                              !
277     !  Purpose: Record information about the location and conditions that  !
278     !           resulted in a negative gas phase volume fraction.          !
279     !                                                                      !
280     !  Author: J. Musser                                  Date: 16-APR-15  !
281     !  Reviewer:                                          Date:            !
282     !                                                                      !
283     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
284           SUBROUTINE EPgErr_LOG(IJK, tHeader)
285     
286     ! Simulation time
287           use run, only: TIME
288     ! Gas phase pressure.
289           use fldvar, only: EP_g, EP_s, ROP_s, RO_s, U_s, V_s, W_s
290           use cutcell
291     
292           use physprop, only: SMAX
293           use compar
294           use indices
295           use functions
296           use error_manager
297     
298     
299           IMPLICIT NONE
300     
301           INTEGER, intent(in) :: IJK
302           LOGICAL, intent(inout) :: tHeader
303     
304           LOGICAL :: lExists
305           CHARACTER(LEN=255) :: lFile
306           INTEGER, parameter :: lUnit = 4868
307           LOGICAL, save :: fHeader = .TRUE.
308     
309           INTEGER :: M
310     
311           lFile = '';
312           if(numPEs > 1) then
313              write(lFile,"('EPgErr_',I4.4,'.log')") myPE
314           else
315              write(lFile,"('EPgErr.log')")
316           endif
317           inquire(file=trim(lFile),exist=lExists)
318           if(lExists) then
319              open(lUnit,file=trim(adjustl(lFile)),                         &
320                 status='old', position='append')
321           else
322              open(lUnit,file=trim(adjustl(lFile)), status='new')
323           endif
324     
325           if(fHeader) then
326              write(lUnit,1000)
327              fHeader = .FALSE.
328           endif
329     
330           if(tHeader) then
331              write(lUnit,"(/2x,'Simulation time: ',g12.5)") TIME
332              tHeader = .FALSE.
333           endif
334     
335           write(lUnit,1001) IJK, I_OF(IJK), J_OF(IJK), K_OF(IJK)
336           write(lUnit,"(6x,A,1X,g12.5)",ADVANCE='NO') 'EP_g:', EP_g(IJK)
337           DO M=1,SMAX
338              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
339                 trim(iVar('EP_s',M)), EP_s(IJK,M)
340           ENDDO
341           write(lUnit,"(' ')")
342     
343           write(lUnit,"(24x)", ADVANCE='NO')
344           DO M=1,SMAX
345              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
346                 trim(iVar('RO_s',M)), RO_s(IJK,M)
347           ENDDO
348           write(lUnit,"(' ')")
349     
350           write(lUnit,"(24x)", ADVANCE='NO')
351           DO M=1,SMAX
352              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
353                 trim(iVar('ROPs',M)), ROP_s(IJK,M)
354           ENDDO
355           write(lUnit,"(24x,' ')")
356     
357     
358           write(lUnit,"(24x)", ADVANCE='NO')
359           DO M=1,SMAX
360              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
361                 trim(iVar('U_se',M)), U_s(IJK,M)
362           ENDDO
363           write(lUnit,"(24x,' ')")
364     
365           write(lUnit,"(24x)", ADVANCE='NO')
366           DO M=1,SMAX
367              write(lUnit,"(2x,A,1X,g12.5)", ADVANCE='NO') &
368                 trim(iVar('V_sn',M)), V_s(IJK,M)
369           ENDDO
370           write(lUnit,"(24x,' ')")
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('W_st',M)), W_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('U_sw',M)), U_s(WEST_OF(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('V_ss',M)), V_s(SOUTH_OF(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('W_sb',M)), W_s(BOTTOM_OF(IJK),M)
397           ENDDO
398           write(lUnit,"(24x,' ')")
399     
400     
401           if(CARTESIAN_GRID) then
402              write(lUnit,"(6x,A,1X,L1)",ADVANCE='NO') 'Cut Cell:', CUT_CELL_AT(IJK)
403              write(lUnit,"(6x,A,1X,L1)") 'Small Cell:', SMALL_CELL_AT(IJK)
404              write(lUnit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
405                 xg_e(I_OF(IJK)), yg_n(J_of(ijk)), zg_t(k_of(ijk))
406           endif
407     
408           close(lUnit)
409     
410           RETURN
411     
412      1000 FORMAT('One or more cells have reported a negative gas volume ', &
413              'fraction (EP_g).',/)
414     
415      1001 FORMAT(/4X,'IJK: ',I8,7X,'I: ',I4,'  J: ',I4,'  K: ',I4)
416     
417           END SUBROUTINE EPgErr_LOG
418