1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv! 2 ! ! 3 ! Subroutine: DRAG_SS_DEM_NONINTERP ! 4 ! ! 5 ! Purpose: This routine is called from the DISCRETE side to calculate ! 6 ! the solids-solids drag force acting on each particle using cell- ! 7 ! center continuum solids velocity. The total contact force is also ! 8 ! updated to include P* for over-packing. ! 9 ! ! 10 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^! 11 SUBROUTINE DRAG_SS_DEM_NONINTERP 12 13 ! Modules 14 !---------------------------------------------------------------------// 15 ! Fluid grid loop bounds. 16 use compar, only: IJKStart3, IJKEnd3 17 ! Fudge factor for SS drag (c.f. Gera, 2004) 18 use constant, only: SEGREGATION_SLOPE_COEFFICIENT 19 ! The count and a list of particles in IJK 20 use discretelement, only: PINC, PIJK 21 use derived_types, only: PIC 22 ! Particle velocity and density 23 use discretelement, only: DES_VEL_NEW, RO_SOL 24 ! Particle radius and volume. 25 use discretelement, only: DES_RADIUS, PVOL 26 ! Total forces acting on particle 27 use discretelement, only: FC 28 ! Gas phase volume fraction 29 use fldvar, only: EP_g 30 ! Solids pressure 31 use fldvar, only: P_STAR 32 ! Diameter of continuum solids phases 33 use fldvar, only: D_P 34 ! Material (ROs) and bulk density (ROPs) of continuum solids 35 use fldvar, only: ROP_s, RO_s 36 use functions, only: IS_NONEXISTENT, IS_GHOST, IS_ENTERING_GHOST, IS_EXITING_GHOST 37 ! Function to deterine if a cell contains fluid. 38 use functions, only: FLUID_AT 39 ! Double precision values. 40 use param1, only: ZERO, ONE 41 ! Array sizes for solids 42 use param, only: DIMENSION_M 43 ! Number of continuum solids phases 44 use physprop, only: SMAX 45 ! Flag that a continuum solids can pack 46 use physprop, only: CLOSE_PACKED 47 ! radial distribution function 48 use rdf, only: g_0 49 IMPLICIT NONE 50 51 ! Local variables 52 !---------------------------------------------------------------------// 53 ! Loop indices: 54 INTEGER :: IJK, M, NINDX, NP, L 55 ! average solids velocity at scalar cell center in array form 56 DOUBLE PRECISION :: VELCS(3, DIMENSION_M) 57 ! relative velocity between solids phase m and l 58 DOUBLE PRECISION :: VSLP(3), VREL 59 ! radial distribution function between phase M and L 60 DOUBLE PRECISION :: G0_ML 61 ! Sum over all phases of ratio volume fraction over particle diameter 62 DOUBLE PRECISION :: EPSoDP 63 ! DEM particle diameter calculated from radius. 64 DOUBLE PRECISION :: lDP 65 ! solid-solid drag coefficient 66 DOUBLE PRECISION :: lDss 67 ! Intermediate calculations for volume fraction. 68 DOUBLE PRECISION :: OoEPg, EPg_2 69 ! Drag force acting on each phase. 70 DOUBLE PRECISION :: D_FORCE(3) 71 72 !......................................................................! 73 74 ! Calculate the solid-solid drag for each particle. 75 !---------------------------------------------------------------------// 76 !!$omp parallel do schedule(guided, 50) default(none) & 77 !!$omp shared(IJKSTART3, IJKEND3, PINC, PIC, DES_VEL_NEW, & 78 !!$omp D_P, RO_s, ROP_s, EP_G, DES_RADIUS, RO_SOL, FC, & 79 !!$omp PVOL, SEGREGATION_SLOPE_COEFFICIENT, CLOSE_PACKED, P_STAR)& 80 !!$omp private(IJK, OoEPg, EPg_2, EPSoDP, NP, lDP, D_FORCE, G0_ML, & 81 !!$omp VELCS, VSLP, VREL, lDss, L) 82 DO IJK = IJKSTART3, IJKEND3 83 84 IF(.NOT.FLUID_AT(IJK)) CYCLE 85 IF(PINC(IJK) == 0) CYCLE 86 87 OoEPg = ONE/EP_g(IJK) 88 EPg_2 = EP_g(IJK)*EP_g(IJK) 89 90 ! Calculate solids volume fraction over particle diameter. 91 CALL CALC_EPSoDP(IJK, EPSoDP) 92 93 ! Calculate the continuum solids velocity at scalar cell center. 94 CALL CALC_CELL_CENTER_CSOLIDS_VEL(IJK, VelCS) 95 96 ! Calculate the solids drag for each particle in the current cell. 97 DO NINDX = 1,PINC(IJK) 98 NP = PIC(IJK)%P(NINDX) 99 ! skipping indices that do not represent particles and ghost particles 100 IF(IS_NONEXISTENT(NP)) CYCLE 101 IF(IS_GHOST(NP) .OR. IS_ENTERING_GHOST(NP) .OR. IS_EXITING_GHOST(NP)) CYCLE 102 103 ! Diameter of particle (not a phase diameter). 104 lDP = 2.0d0*DES_RADIUS(NP) 105 106 D_FORCE = ZERO 107 DO M = 1, SMAX 108 109 ! evaluating g0 - taken from G_0.f subroutine (lebowitz form) 110 ! this section is needed to account for all solids phases until g0 for 111 ! multiple solids types (i.e. discrete & continuum) can be addressed 112 ! more effectively. 113 ! G0_ML = OoEPg + 3.0d0*EPSoDP*D_P(IJK,M)*lDP / & 114 ! (EPg_2 *(D_P(IJK,M) + lDP)) 115 ! the calculation below allows for a more generic g0 form, however, 116 ! the calculation above is truer to the discrete nature of the info 117 L = PIJK(NP,5) 118 G0_ML = G_0(IJK,L,M) 119 120 ! Relative (slip) velocity. 121 VSLP = DES_VEL_NEW(NP,:) - VelCS(:,M) 122 ! Relative velocity magnitude. 123 VREL = sqrt(dot_product(VSLP,VSLP)) 124 125 CALL DRAG_SS_SYAM0(lDss, D_P(IJK,M), lDP, RO_S(IJK,M), & 126 RO_SOL(NP), G0_ML, VREL) 127 128 lDss = lDss*ROP_S(IJK,M)*RO_Sol(NP) 129 130 ! accounting for particle-particle drag due to enduring contact in a 131 ! close-packed system 132 IF(CLOSE_PACKED(M)) lDss = lDss + & 133 SEGREGATION_SLOPE_COEFFICIENT*P_star(IJK) 134 135 ! Calculating the accumulated solids-solids drag force. 136 D_FORCE(:) = D_FORCE(:) - lDss*VSLP(:) 137 138 ENDDO ! end do loop (M=1,SMAX) 139 140 FC(NP,:) = FC(NP,:) + D_FORCE(:)*PVOL(NP) 141 142 ENDDO ! END DO LOOP (NP=1,MAX_PIP) 143 ENDDO ! END DO LOOP (IJK=IJKSTART3, IJKEND3) 144 !!$omp end parallel do 145 146 RETURN 147 END SUBROUTINE DRAG_SS_DEM_NONINTERP 148 149 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv! 150 ! ! 151 ! Subroutine: DES_DRAG_SS ! 152 ! ! 153 ! Purpose: This subroutine is called from the CONTINUUM side. It ! 154 ! calculate the solids-solids drag force coefficient between ! 155 ! continuum and discrete solids. ! 156 ! ! 157 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^! 158 SUBROUTINE DRAG_SS_TFM_NONINTERP 159 160 ! Modules 161 !---------------------------------------------------------------------// 162 ! Fluid grid loop bounds. 163 use compar, only: IJKStart3, IJKEnd3 164 ! Fudge factor for SS drag (c.f. Gera, 2004) 165 use constant, only: SEGREGATION_SLOPE_COEFFICIENT 166 ! The count and a list of particles in IJK 167 use derived_types, only: PIC 168 use discretelement, only: PINC, PIJK, DES_VOL_NODE 169 ! Particle velocity and density 170 use discretelement, only: DES_VEL_NEW, RO_SOL 171 ! Particle radius and volume. 172 use discretelement, only: DES_RADIUS 173 ! Total forces acting on particle 174 use discretelement, only: SDRAG_AM, SDRAG_BM, F_SDS 175 ! Gas phase volume fraction 176 use fldvar, only: EP_g 177 ! Diameter of continuum solids phases 178 use fldvar, only: D_P 179 ! Material (ROs) and bulk density (ROPs) of continuum solids 180 use fldvar, only: ROP_s, RO_s 181 ! Solids pressure 182 use fldvar, only: P_STAR 183 ! Function to deterine if a cell contains fluid. 184 use functions, only: FLUID_AT 185 use functions, only: IS_NONEXISTENT, IS_GHOST 186 use functions, only: IS_ENTERING_GHOST, IS_EXITING_GHOST 187 use geometry, only: vol 188 ! Double precision values. 189 use param1, only: ZERO, ONE 190 ! Array sizes for solids 191 use param, only: DIMENSION_M 192 ! Flag that a continuum solids can pack 193 use physprop, only: CLOSE_PACKED 194 ! Number of continuum solids phases 195 use physprop, only: SMAX 196 ! radial distribution function 197 use rdf, only: g_0 198 IMPLICIT NONE 199 200 ! Local variables 201 !---------------------------------------------------------------------// 202 ! Loop indices: 203 INTEGER :: IJK, M, NINDX, NP, L 204 ! average solids velocity at scalar cell center in array form 205 DOUBLE PRECISION :: VELCS(3, DIMENSION_M) 206 ! relative velocity between solids phase m and l 207 DOUBLE PRECISION :: VSLP(3), VREL 208 ! radial distribution function between phase M and L 209 DOUBLE PRECISION :: G0_ML 210 ! Sum over all phases of ratio volume fraction over particle diameter 211 DOUBLE PRECISION :: EPSoDP 212 ! DEM particle diameter calculated from radius. 213 DOUBLE PRECISION :: lDP 214 ! solid-solid drag coefficient 215 DOUBLE PRECISION :: lDss 216 ! Intermediate calculations for volume fraction. 217 DOUBLE PRECISION :: OoEPg, EPg_2 218 ! Drag force acting on each phase. 219 DOUBLE PRECISION :: lFORCE 220 !......................................................................! 221 222 DO IJK = IJKSTART3, IJKEND3 223 224 F_SDS(IJK,:) = ZERO 225 SDRAG_AM(IJK,:) = ZERO 226 SDRAG_BM(IJK,:,:) = ZERO 227 228 IF(.NOT.FLUID_AT(IJK)) CYCLE 229 IF(PINC(IJK) == 0) CYCLE 230 231 OoEPg = ONE/EP_g(IJK) 232 EPg_2 = EP_g(IJK)*EP_g(IJK) 233 234 ! Calculate solids volume fraction over particle diameter. 235 CALL CALC_EPSoDP(IJK, EPSoDP) 236 237 ! Calculate the continuum solids velocity at scalar cell center. 238 CALL CALC_CELL_CENTER_CSOLIDS_VEL(IJK, VelCS) 239 240 ! Calculate the solids drag for each particle in the current cell. 241 DO NINDX = 1,PINC(IJK) 242 NP = PIC(IJK)%P(NINDX) 243 ! skipping indices that do not represent particles and ghost particles 244 IF(IS_NONEXISTENT(NP)) CYCLE 245 IF(IS_GHOST(NP) .OR. IS_ENTERING_GHOST(NP) .OR. IS_EXITING_GHOST(NP)) CYCLE 246 247 ! Diameter of particle (not a phase diameter). 248 lDP = 2.0d0*DES_RADIUS(NP) 249 250 DO M = 1, SMAX 251 252 ! Relative (slip) velocity. 253 VSLP = DES_VEL_NEW(NP,:) - VelCS(:,M) 254 ! Relative velocity magnitude. 255 VREL = sqrt(dot_product(VSLP,VSLP)) 256 257 ! evaluating g0 - taken from G_0.f subroutine (lebowitz form) 258 ! this section is needed to account for all solids phases until g0 for 259 ! multiple solids types (i.e. discrete & continuum) can be addressed 260 ! more effectively. 261 ! G0_ML = OoEPg + 3.0d0*EPSoDP*D_P(IJK,M)*lDP / & 262 ! (EPg_2 *(D_P(IJK,M) + lDP)) 263 ! the calculation below allows for a more generic g0 form, however, 264 ! the calculation above is truer to the discrete nature of the info 265 L = PIJK(NP,5) 266 G0_ML = G_0(IJK,L,M) 267 268 CALL DRAG_SS_SYAM0(lDss, D_P(IJK,M), lDP, RO_S(IJK,M), & 269 RO_SOL(NP), G0_ML, VREL) 270 271 lDss = lDss*ROP_S(IJK,M)*RO_Sol(NP) 272 273 ! accounting for particle-particle drag due to enduring contact in a 274 ! close-packed system 275 IF(CLOSE_PACKED(M)) lDss = lDss + & 276 SEGREGATION_SLOPE_COEFFICIENT*P_star(IJK) 277 278 ! Calculating the accumulated solids-solids drag force. 279 lFORCE = lDss/vol(L) 280 281 SDRAG_AM(IJK,M) = SDRAG_AM(IJK,M) + lFORCE 282 SDRAG_BM(IJK,:,M) = SDRAG_BM(IJK,:,M) + & 283 lFORCE*DES_VEL_NEW(NP,:) 284 285 ENDDO ! end do loop (M=1,SMAX) 286 287 ENDDO ! END DO LOOP (NP=1,MAX_PIP) 288 289 F_SDS(IJK,:) = SDRAG_AM(IJK,:) 290 291 ENDDO ! END DO LOOP (IJK=IJKSTART3, IJKEND3) 292 293 RETURN 294 END SUBROUTINE DRAG_SS_TFM_NONINTERP 295 296 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv! 297 ! ! 298 ! Subroutine: CALC_CELL_CENTER_CSOLIDS_VEL ! 299 ! Author: J.Musser Date: 07-NOV-14 ! 300 ! ! 301 ! Purpose: Calculate the scalar cell center continuum solids ! 302 ! velocity. This code is common to the DEM and GAS calls for non- ! 303 ! interpolated drag routines. ! 304 ! ! 305 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^! 306 SUBROUTINE CALC_CELL_CENTER_CSOLIDS_VEL(IJK, lVELFP) 307 308 ! Modules 309 !---------------------------------------------------------------------// 310 ! Flags and correction factors for cut momentum cells. 311 use cutcell, only: CUT_U_TREATMENT_AT, THETA_UE, THETA_UE_BAR 312 use cutcell, only: CUT_V_TREATMENT_AT, THETA_VN, THETA_VN_BAR 313 use cutcell, only: CUT_W_TREATMENT_AT, THETA_WT, THETA_WT_BAR 314 ! Solids velocities. 315 use fldvar, only: U_S, V_S, W_S 316 ! Functions to average momentum to scalar cell center. 317 use fun_avg, only: AVG_X_E, AVG_Y_N, AVG_Z_T 318 ! Functions to lookup adjacent cells by index. 319 use functions, only: IM_OF, JM_OF, KM_OF 320 ! Flag for 3D simulatoins. 321 use geometry, only: DO_K 322 use indices, only: I_OF 323 ! Double precision parameters. 324 use param1, only: ZERO 325 ! Array sizes for solids 326 use param, only: DIMENSION_M 327 ! Number of solid phases. 328 use physprop, only: SMAX 329 IMPLICIT NONE 330 331 ! Dummy arguments 332 !---------------------------------------------------------------------// 333 ! Fluid cell and solids phase indices 334 INTEGER, INTENT(IN) :: IJK 335 ! Fluid velocity vector at IJK cell center. 336 DOUBLE PRECISION, INTENT(OUT) :: lVELFP(3, DIMENSION_M) 337 338 ! Local variables 339 !---------------------------------------------------------------------// 340 ! Phase loop counter 341 INTEGER :: M 342 ! Indices of adjacent cells 343 INTEGER :: IMJK, IJMK, IJKM 344 !......................................................................! 345 346 ! Calculate the average fluid velocity at scalar cell center. 347 DO M=1,SMAX 348 349 IMJK = IM_OF(IJK) 350 IF(CUT_U_TREATMENT_AT(IMJK)) THEN 351 lVELFP(1,M) = (THETA_UE_BAR(IMJK)*U_S(IMJK,M) + & 352 THETA_UE(IMJK)*U_S(IJK,M)) 353 ELSE 354 lVELFP(1,M) = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I_OF(IJK)) 355 ENDIF 356 357 IJMK = JM_OF(IJK) 358 IF(CUT_V_TREATMENT_AT(IJMK)) THEN 359 lVELFP(2,M) = (THETA_VN_BAR(IJMK)*V_S(IJMK,M) + & 360 THETA_VN(IJMK)*V_S(IJK,M)) 361 ELSE 362 lVELFP(2,M) = AVG_Y_N(V_S(IJMK,M),V_S(IJK,M)) 363 ENDIF 364 365 IF(DO_K) THEN 366 IJKM = KM_OF(IJK) 367 IF(CUT_W_TREATMENT_AT(IJKM)) THEN 368 lVELFP(3,M) = (THETA_WT_BAR(IJKM)*W_S(IJKM,M) + & 369 THETA_WT(IJKM)* W_S(IJK,M)) 370 ELSE 371 lVELFP(3,M) = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M)) 372 ENDIF 373 ELSE 374 lVELFP(3,M) = ZERO 375 ENDIF 376 ENDDO 377 378 RETURN 379 END SUBROUTINE CALC_CELL_CENTER_CSOLIDS_VEL 380 381 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv! 382 ! ! 383 ! Subroutine: CALC_EPSoDP ! 384 ! Author: J.Musser Date: 07-NOV-14 ! 385 ! ! 386 ! ! 387 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^! 388 SUBROUTINE CALC_EPSoDP(IJK, lEPSoDP) 389 390 ! Global Variables: 391 !---------------------------------------------------------------------// 392 ! The count and a list of particles in IJK 393 use derived_types, only: PIC 394 use discretelement, only: PINC 395 ! Particle volume and radius 396 use discretelement, only: PVOL, DES_RADIUS 397 ! Volume of scalar cell 398 use geometry, only: VOL 399 ! Continuum solids diamter 400 use fldvar, only: D_P 401 ! Function to calculate continuum solids volume fraction 402 use fldvar, only: EP_s 403 ! Number of continuum solids phases 404 use physprop, only: SMAX 405 406 use functions 407 408 ! Global Parameters: 409 !---------------------------------------------------------------------// 410 ! Double precision parameters. 411 use param1, only: ZERO 412 413 IMPLICIT NONE 414 415 ! Dummy arguments 416 !---------------------------------------------------------------------// 417 ! Fluid cell and solids phase indices 418 INTEGER, INTENT(IN) :: IJK 419 ! Fluid velocity vector at IJK cell center. 420 DOUBLE PRECISION, INTENT(OUT) :: lEPSoDP 421 422 ! Local variables: 423 !---------------------------------------------------------------------// 424 ! Loop counters 425 INTEGER :: NP, NINDX, M 426 427 ! Calculate the sum of the particle volumes divided by radius. 428 lEPSoDP = ZERO 429 DO NINDX = 1,PINC(IJK) 430 NP = PIC(IJK)%P(NINDX) 431 IF(IS_NONEXISTENT(NP)) CYCLE 432 IF(IS_GHOST(NP) .OR. IS_ENTERING_GHOST(NP) .OR.& 433 IS_EXITING_GHOST(NP)) CYCLE 434 lEPSoDP = lEPSoDP + PVOL(NP)/DES_RADIUS(NP) 435 ENDDO 436 ! Convert radius to diameter and divide by cell volume. 437 lEPSoDP = lEPSoDP/(2.0d0*VOL(IJK)) 438 439 ! Add contributions from continuum solids. 440 DO M = 1, SMAX 441 lEPSoDP = lEPSoDP + EP_s(IJK,M)/D_p(IJK,M) 442 ENDDO 443 444 RETURN 445 END SUBROUTINE CALC_EPSoDP 446