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