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