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 ! Reviewer: Date: C 16 ! C 17 ! C 18 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 19 20 SUBROUTINE CALC_VOL_FR(P_STAR, RO_G, ROP_G, EP_G, ROP_S, IER) 21 22 !----------------------------------------------- 23 ! Modules 24 !----------------------------------------------- 25 USE param 26 USE param1 27 USE parallel 28 USE run 29 USE geometry 30 USE indices 31 USE physprop 32 USE visc_s 33 USE constant 34 USE pgcor 35 USE pscor 36 USE compar 37 USE sendrecv 38 USE discretelement 39 USE mpi_utility 40 use fldvar, only: RO_S, EP_S 41 USE solids_pressure 42 USE functions 43 44 IMPLICIT NONE 45 !----------------------------------------------- 46 ! Dummy Arguments 47 !----------------------------------------------- 48 ! Solids pressure 49 DOUBLE PRECISION, INTENT(IN) :: P_star(DIMENSION_3) 50 ! Gas density 51 DOUBLE PRECISION, INTENT(INOUT) :: RO_g(DIMENSION_3) 52 ! Gas bulk density 53 DOUBLE PRECISION, INTENT(INOUT) :: ROP_g(DIMENSION_3) 54 ! Gas volume fraction 55 DOUBLE PRECISION, INTENT(INOUT) :: EP_g(DIMENSION_3) 56 ! solids bulk densities 57 DOUBLE PRECISION, INTENT(INOUT) :: ROP_s(DIMENSION_3, DIMENSION_M) 58 ! error index 59 INTEGER, INTENT(INOUT) :: IER 60 !----------------------------------------------- 61 ! Local variables 62 !----------------------------------------------- 63 ! volume of particle type M for GHD theory 64 DOUBLE PRECISION :: VOL_M 65 ! volume fraction of close-packed region 66 DOUBLE PRECISION :: EPcp 67 ! volume fraction of solids phase 68 DOUBLE PRECISION :: EPS 69 ! sum of volume fractions 70 DOUBLE PRECISION :: SUMVF 71 ! Whichever phase is given by MF becomes the phase whose volume fraction 72 ! is corrected based on all other phases present. Generally MF gets 73 ! defined as 0 so that the gas phase void fraction becomes corrected 74 ! based on the summation of the volume fractions of all other phases 75 INTEGER :: MF 76 ! Whichever phase is given by MCPl becomes the phase whose volume 77 ! fraction is corrected based on value of maximum close packing and all 78 ! other solids phase that can close pack. This is basically a local 79 ! variable for MCP but only in cells where close packed conditions 80 ! exist. 81 INTEGER :: MCPl 82 ! Index of solids phase 83 INTEGER :: M 84 ! Indices 85 INTEGER :: IJK 86 87 ! Arrays for storing errors: 88 ! 110 - Negative volume fraciton 89 ! 11x - Unclassified 90 INTEGER :: Err_l(0:numPEs-1) ! local 91 INTEGER :: Err_g(0:numPEs-1) ! global 92 93 !----------------------------------------------- 94 95 ! Initialize error flag. 96 Err_l = 0 97 98 !!$omp parallel do private(MCPl, EPCP, SUMVF, MF, M) & 99 !!$omp& schedule(static) 100 101 DO IJK = ijkstart3, ijkend3 102 IF (FLUID_AT(IJK)) THEN 103 104 ! calculate the volume fraction of the solids phase based on the volume 105 ! fractions of all other solids phases and on the value of maximum 106 ! packing when that solids phase continuity was not solved for the 107 ! indicated cell. 108 !----------------------------------------------------------------->>> 109 IF (PHASE_4_P_S(IJK) /= UNDEFINED_I) THEN 110 ! for the current cell check if a solids phase continuity was skipped. 111 ! this value will either be undefined (no cell was skipped in any 112 ! solids continuity equations) or be a solids phase index (indicated 113 ! solids continuity was skipped) for a solids phase that does close pack 114 ! (i.e., close_packed=T) and the cell exhibits close packing conditions. 115 ! note: this branch is never entered given the existing version of 116 ! mark_phase_4_cor. 117 MCPl = PHASE_4_P_s(IJK) 118 119 ! finding the value of maximum close packing based on the expression for 120 ! plastic pressure 121 EPCP = 1. - INV_H(P_STAR(IJK),EP_g_blend_end(ijk)) 122 ! summing the solids volume fraction of all continuum solids phases 123 ! that are marked as close_packed except for the solids phase which 124 ! is also marked by phase_4_p_s (skipping its own continuity) 125 SUMVF = ZERO 126 DO M = 1, MMAX 127 IF (CLOSE_PACKED(M) .AND. M/=MCPl) SUMVF = SUMVF + EP_S(IJK,M) 128 ENDDO 129 ! sum of solids volume fraction from the DEM 130 IF (DES_CONTINUUM_HYBRID) THEN 131 DO M = 1,DES_MMAX 132 EPS = DES_ROP_S(IJK,M)/DES_RO_S(M) 133 SUMVF = SUMVF + EPS 134 ENDDO 135 ENDIF 136 ROP_S(IJK,MCPl) = (EPCP - SUMVF)*RO_S(IJK,MCPl) 137 ENDIF 138 !-----------------------------------------------------------------<<< 139 140 141 ! calculate the volume fraction of the 'solids' phase based on the 142 ! volume fractions of all other phases (including gas) if the gas 143 ! continuity was solved rather than that phases own continuity. 144 ! if the gas continuity was not solved then calculate the void 145 ! fraction of the gas phase based on all other solids phases. 146 !----------------------------------------------------------------->>> 147 ! for the current cell check if the gas phase continuity was solved for 148 ! the gas phase while the solids phase continuity for the indicated 149 ! solids phase was skipped. this value will either be 0 (gas phase 150 ! continuity was not solved) or be a solids phase index (gas continuity 151 ! was solved while the indicated solids phase continuity was skipped) 152 ! for a solids phase that does not close pack (i.e., close_packed=F) and 153 ! is in greater concentration than the gas phase. 154 ! Note: MF will always be set to 0 here given the existing version of 155 ! mark_phase_4_cor 156 MF = PHASE_4_P_G(IJK) 157 158 SUMVF = ZERO 159 IF (0 /= MF) THEN 160 ! if gas continuity was solved rather than the solids phase, then 161 ! include the gas phase void fraction in the summation here. 162 EP_G(IJK) = ROP_G(IJK)/RO_G(IJK) 163 SUMVF = SUMVF + EP_G(IJK) 164 ENDIF 165 166 ! modified for GHD theory 167 IF(KT_TYPE_ENUM == GHD_2007) THEN 168 ROP_S(IJK,MMAX) = ZERO ! mixture density 169 DO M = 1, SMAX 170 ! volume of particle M based on fixed diamter Dp0 171 VOL_M = PI*D_P0(M)**3/6d0 172 IF (M /= MF) THEN 173 SUMVF = SUMVF + EP_S(IJK,M) 174 ROP_S(IJK,MMAX) = ROP_S(IJK,MMAX) + RO_S(IJK,M)*EP_S(IJK,M) 175 ENDIF 176 ENDDO 177 ELSE 178 ! summing the solids volume fraction of all continuum solids phases 179 ! except for the continuum solids phase which was marked by phase_4_p_g 180 ! (skipping its continuity while solving gas continuity) 181 DO M = 1, MMAX 182 IF (M /= MF) SUMVF = SUMVF + EP_S(IJK,M) 183 ENDDO 184 ENDIF ! end if/else trim(kt_type)=='ghd' 185 186 ! sum of solids volume fraction from the DEM 187 IF (DES_CONTINUUM_HYBRID) THEN 188 DO M = 1,DES_MMAX 189 EPS = DES_ROP_S(IJK,M)/DES_RO_S(M) 190 SUMVF = SUMVF + EPS 191 ENDDO 192 ENDIF 193 194 IF (0 == MF) THEN 195 ! if no gas phase continuity was solved in the current cell then correct 196 ! the void fraction of the gas phase based on the total solids volume 197 ! fraction of all solids phases 198 EP_G(IJK) = ONE - SUMVF 199 200 ! Set error flag for negative volume fraction. 201 IF (EP_G(IJK) < ZERO) Err_l(myPE) = 110 202 203 ROP_G(IJK) = EP_G(IJK)*RO_G(IJK) 204 ELSE 205 ! else correct the volume fraction of the solids phase that was marked 206 ROP_S(IJK,MF) = (ONE - SUMVF)*RO_S(IJK,MF) 207 ENDIF 208 !-----------------------------------------------------------------<<< 209 210 ENDIF ! end if (fluid_at(ijk)) 211 ENDDO ! end do (ijk=ijkstart3,ijkend3) 212 213 CALL send_recv(EP_G, 2) 214 CALL send_recv(ROP_G, 2) 215 CALL send_recv(ROP_S, 2) 216 217 CALL global_all_sum(Err_l, Err_g) 218 IER = maxval(Err_g) 219 220 221 RETURN 222 END SUBROUTINE CALC_VOL_FR 223 224 225