1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 2 ! C 3 ! Module name: DISCRETELEMENT C 4 ! Purpose: DES mod file C 5 ! Common Block containing DEM conditions C 6 ! C 7 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 8 9 MODULE DISCRETELEMENT 10 11 !----------------------------------------------- 12 ! Modules 13 !----------------------------------------------- 14 USE param 15 USE param1 16 IMPLICIT NONE 17 !----------------------------------------------- 18 19 ! Define interface - needed when passing arrays of assumed size 20 ! This function is used to identify fluid grid index (i,j or k) of a 21 ! particle in a given x-, y-, z- line of space 22 INTERFACE 23 INTEGER FUNCTION DES_GETINDEXFROMPOS(LIM1,LIM2,PART_POS,& 24 GRID_POS,AXIS,AXIS_INDEX) 25 INTEGER :: LIM1, LIM2 26 DOUBLE PRECISION :: PART_POS 27 DOUBLE PRECISION, DIMENSION(:) :: GRID_POS 28 CHARACTER(LEN=1) :: AXIS,AXIS_INDEX 29 END FUNCTION DES_GETINDEXFROMPOS 30 END INTERFACE 31 32 33 ! Total number of particles in simulation: read from input or generated 34 INTEGER PARTICLES 35 36 ! Start particle tracking quantities 37 !----------------------------------------------------------------->>> 38 ! Generally for inlet/outlet related routines but also employed in 39 ! tracking for parallelization 40 41 ! Dynamic particle states: 42 43 ! NONEXISTENT: This state is used with the inlet/outlet to skip 44 ! indices that do not represent particles in the system or indices 45 ! that represent particles that have exited the system. 46 47 ! ENTERING: This state identifies a particle as 'new' if true. 48 ! Particles with a classification of 'new' do not react when in contact 49 ! with a wall or another particle, however existing particles do collide 50 ! and interact with 'new' particles. The classification allows new 51 ! particles to push particles already in the system out of the way when 52 ! entering to prevent overlap. This flag is also used when the center 53 ! of a particle crosses a dem outlet (i.e. an exiting particle; see 54 ! EXITING) so that the particle will maintain its present trajectory 55 ! until it has fully exited the system 56 57 ! EXITING: This state identifies a particle as 'exiting' if true. 58 ! If a particle initiates contact with a wall surface designated as a 59 ! des outlet, this flag is set to true. With this classification the 60 ! location of the particle is checked to assess if the particle has 61 ! fully exited the system. At this point, the particle is removed 62 ! from the list. 63 64 ! GHOST, ENTERING_GHOST, EXITING_GHOST: for ghost particles 65 66 INTEGER(KIND=1), DIMENSION(:), ALLOCATABLE :: PARTICLE_STATE ! (PARTICLES) 67 68 INTEGER, PARAMETER :: nonexistent=0 69 INTEGER, PARAMETER :: normal_particle=1 70 INTEGER, PARAMETER :: entering_particle=2 71 INTEGER, PARAMETER :: exiting_particle=3 72 INTEGER, PARAMETER :: normal_ghost=4 73 INTEGER, PARAMETER :: entering_ghost=5 74 INTEGER, PARAMETER :: exiting_ghost=6 75 76 ! PARALLEL PROCESSING: explanation of variables in parallel architecture 77 ! pip - particles in each processor (includes the ghost particles) 78 ! max_pip - maximum allocated particles in processor 79 80 ! Number of particles in the system (current) 81 INTEGER PIP 82 ! Global sum of particles (excluding ghost) in the system 83 INTEGER TOT_PAR 84 ! Maximum particles permitted in the system at once 85 INTEGER MAX_PIP 86 87 ! End particle tracking quantities 88 !-----------------------------------------------------------------<<< 89 90 ! For parallel processing: global id of particles 91 INTEGER, DIMENSION(:), ALLOCATABLE :: IGLOBAL_ID 92 ! Ghost count of particles on each processor 93 INTEGER :: IGHOST_CNT 94 ! Maximum global id, new particles global id will be assigned based 95 ! on this value 96 Integer :: imax_global_id 97 98 99 ! If gener_part_config is true, then the particle_input.dat file does 100 ! not need to be supplied nor does the total number of particles as 101 ! these are determined based on the specified volume fraction (vol_frac) 102 ! in the specified domain (des_eps_xyzstart) 103 LOGICAL :: GENER_PART_CONFIG 104 DOUBLE PRECISION :: VOL_FRAC(DIM_M) 105 DOUBLE PRECISION :: DES_EPS_XSTART, & 106 DES_EPS_YSTART, DES_EPS_ZSTART 107 ! volume of the IC region for computing number of particles to be seeded 108 DOUBLE PRECISION, dimension(:), allocatable :: VOL_IC_REGION!(DIMENSION_IC) 109 ! number of particles for each phase corresponding to the IC number. This will 110 ! be real particles for DEM but parcels or computational particles for PIC model 111 INTEGER, dimension(:,:), allocatable :: PART_MPHASE_BYIC!(DIMENSION_IC, DIM_M) 112 113 ! Number of real particles by IC and by solid phase. Only relevant for PIC model 114 double precision, dimension(:,:), allocatable :: REALPART_MPHASE_BYIC!(DIMENSION_IC, DIM_M) 115 ! The number of particles that belong to solid phase M according to the 116 ! vol_frac and particle diameter. this information is used when 117 ! gener_part_config is invoked for initialization 118 ! This will be removed soon as PART_MPHASE_BYIC will be used from now on 119 INTEGER PART_MPHASE(DIM_M) 120 121 ! The number of real particles that belong to solid phase M during the initialization. 122 ! It is equal to Part_mphase for DEM but implies real number of particles for PIC model 123 double precision REALPART_MPHASE(DIM_M) 124 ! Assigns the initial particle velocity distribution based on user 125 ! specified mean and standard deviation (regardless if already set 126 ! within particle_input.dat) 127 DOUBLE PRECISION pvel_mean, PVEL_StDev 128 129 ! Output/debug controls 130 !----------------------------------------------------------------->>> 131 ! Logic that controls whether to print data dem simulations (granular or 132 ! coupled) 133 LOGICAL PRINT_DES_DATA 134 135 ! logic that controls if des run time messages are printed on screen or not 136 LOGICAL PRINT_DES_SCREEN 137 138 ! This specifies the file type used for outputting DES data 139 ! options are : 140 ! TECPLOT - data is written in Tecplot format 141 ! undefined - data is written in ParaView format (default) 142 CHARACTER(LEN=64) :: DES_OUTPUT_TYPE 143 144 ! Used sporadically to control screen dumps (for debug purposes) 145 LOGICAL :: DEBUG_DES 146 147 ! Single particle no. index that is followed if debugging 148 INTEGER FOCUS_PARTICLE 149 150 ! Output file count for .vtp type files and for tecplot files; 151 ! for vtp output used to label .vtp file names in sequential order 152 ! and is saved so restarts begin at the correct count 153 INTEGER VTP_FINDEX, TECPLOT_FINDEX 154 ! End Output/debug controls 155 !-----------------------------------------------------------------<<< 156 157 ! DES - Continuum 158 LOGICAL DISCRETE_ELEMENT 159 LOGICAL DES_CONTINUUM_COUPLED 160 161 ! DES - Invoke hybrid model where both the DEM and continuum model 162 ! are employed to describe solids 163 LOGICAL DES_CONTINUUM_HYBRID 164 165 ! DES - 166 ! With this logic the particles see the fluid but the fluid does 167 ! not see the particles. 168 LOGICAL DES_ONEWAY_COUPLED 169 170 ! Only used when coupled and represents the number of times a pure 171 ! granular flow simulation is run before the actual coupled simulation 172 ! is started (i.e. for particle settling w/o fluid forces) 173 INTEGER NFACTOR 174 175 ! Drag 176 LOGICAL TSUJI_DRAG 177 178 ! Collision model, options are as follows 179 ! linear spring dashpot model (default/undefined) 180 ! 'hertzian' model 181 CHARACTER(LEN=64) :: DES_COLL_MODEL 182 INTEGER DES_COLL_MODEL_ENUM 183 INTEGER,PARAMETER :: HERTZIAN=0 184 INTEGER,PARAMETER :: LSD=1 185 186 ! Integration method, options are as follows 187 ! 'euler' first-order scheme (default) 188 ! 'adams_bashforth' second-order scheme (by T.Li) 189 CHARACTER(LEN=64) :: DES_INTG_METHOD 190 LOGICAL INTG_ADAMS_BASHFORTH 191 LOGICAL INTG_EULER 192 193 ! Value of solids time step based on particle properties 194 DOUBLE PRECISION DTSOLID 195 ! Run time value of simulation time used in dem simulation 196 DOUBLE PRECISION S_TIME 197 198 199 ! Neighbor search related quantities 200 !----------------------------------------------------------------->>> 201 ! Neighbor search method, options are as follows 202 ! 1= nsquare, 2=quadtree, 3=octree, 4=grid/cell based search 203 INTEGER DES_NEIGHBOR_SEARCH 204 205 ! Quantities used to determine whether neighbor search should be called 206 INTEGER NEIGHBOR_SEARCH_N 207 DOUBLE PRECISION NEIGHBOR_SEARCH_RAD_RATIO 208 LOGICAL DO_NSEARCH 209 210 ! Flag on whether to have DES_*_OLD arrays, if either Adams Bashforth or PIC is used 211 LOGICAL DO_OLD 212 213 ! Factor muliplied by sum of radii in grid based neighbor search and 214 ! nsquare search method. increases the effective radius of a particle 215 ! for detecting particle contacts 216 DOUBLE PRECISION FACTOR_RLM 217 218 ! Stores number of neighbors based on neighbor search 219 INTEGER, DIMENSION(:), ALLOCATABLE :: NEIGHBOR_INDEX 220 INTEGER, DIMENSION(:), ALLOCATABLE :: NEIGHBOR_INDEX_OLD 221 INTEGER, DIMENSION(:), ALLOCATABLE :: NEIGHBORS 222 INTEGER, DIMENSION(:), ALLOCATABLE :: NEIGHBORS_OLD 223 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: PFT_NEIGHBOR 224 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: PFT_NEIGHBOR_OLD 225 226 INTEGER, DIMENSION(:), ALLOCATABLE :: CELLNEIGHBOR_FACET_NUM, CELLNEIGHBOR_FACET_MAX 227 INTEGER :: NEIGH_NUM,NEIGH_MAX 228 229 ! Quantities used for reporting: max no. neighbors and max overlap 230 ! that exists during last solid time step of dem simulation 231 DOUBLE PRECISION OVERLAP_MAX 232 233 ! The number of i, j, k divisions in the grid used to perform the 234 ! cell based neighbor search 235 INTEGER :: DESGRIDSEARCH_IMAX, DESGRIDSEARCH_JMAX, & 236 DESGRIDSEARCH_KMAX 237 238 ! End neighbor search related quantities 239 !-----------------------------------------------------------------<<< 240 241 242 ! User specified dimension of the system (by default 2D, but if 3D system is 243 ! desired then it must be explicitly specified) 244 INTEGER, PARAMETER :: DIMN = 3 245 246 ! Variable that is set to the number of walls in the system 247 INTEGER NWALLS 248 249 ! Position of domain boundaries generally given as 250 ! (0, xlength, 0, ylength, 0, zlength) 251 DOUBLE PRECISION WX1, EX2, BY1, TY2, SZ1, NZ2 252 253 ! X, Y, Z position of cell faces of computational fluid grid 254 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: XE !(0:DIMENSION_I) 255 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: YN !(0:DIMENSION_J) 256 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: ZT !(0:DIMENSION_K) 257 258 ! Wall normal vector 259 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: WALL_NORMAL !(NWALLS,3) 260 261 ! Gravity vector and magnitude 262 DOUBLE PRECISION :: GRAV(3) 263 DOUBLE PRECISION :: GRAV_MAG 264 265 266 ! Periodic wall BC 267 LOGICAL DES_PERIODIC_WALLS 268 LOGICAL DES_PERIODIC_WALLS_X 269 LOGICAL DES_PERIODIC_WALLS_Y 270 LOGICAL DES_PERIODIC_WALLS_Z 271 272 273 ! Lees & Edwards wall BC (lost in current DEM) 274 !----------------------------------------------------------------->>> 275 ! Logic for Lees & Edwards BC (T = turn on LE BC) 276 LOGICAL DES_LE_BC 277 ! Relative velocity of LE boundaries (distance/time) 278 DOUBLE PRECISION DES_LE_REL_VEL 279 ! Shear direction 280 ! 2D options are DUDY or DVDX 281 ! 3D options are DUDY, DUDZ, DVDX, DVDZ, DWDX or DWDY 282 ! Note that all other directions are treated as periodic boundaries 283 CHARACTER(LEN=4) :: DES_LE_SHEAR_DIR 284 ! End LE BC 285 !-----------------------------------------------------------------<<< 286 287 288 289 ! Particle-particle and Particle-wall collision model parameters 290 !----------------------------------------------------------------->>> 291 ! Spring contants 292 DOUBLE PRECISION KN, KN_W !Normal 293 DOUBLE PRECISION KT, KT_W, KT_FAC, KT_W_FAC ! Tangential factors = KT/KN and KT_w/KN_w, resp. 294 295 ! Damping coeffients 296 DOUBLE PRECISION ETA_DES_N, ETA_N_W !Normal 297 DOUBLE PRECISION ETA_DES_T, ETA_T_W !Tangential 298 299 ! Flag to use van der Hoef et al. (2006) model for adjusting the rotation of the 300 ! contact plane 301 LOGICAL :: USE_VDH_DEM_MODEL 302 ! Tangential damping factors, eta_t = eta_t_factor * eta_N 303 DOUBLE PRECISION DES_ETAT_FAC, DES_ETAT_W_FAC 304 305 ! Damping coeffients in array form 306 DOUBLE PRECISION :: DES_ETAN(DIM_M, DIM_M), DES_ETAN_WALL(DIM_M) 307 DOUBLE PRECISION :: DES_ETAT(DIM_M, DIM_M), DES_ETAT_WALL(DIM_M) 308 309 ! Friction coeficients 310 DOUBLE PRECISION MEW, MEW_W 311 312 ! coeff of restituion input in one D array, solid solid 313 ! Tangential rest. coef. are used for hertzian collision model but not linear 314 DOUBLE PRECISION DES_EN_INPUT(DIM_M+DIM_M*(DIM_M-1)/2) 315 DOUBLE PRECISION DES_ET_INPUT(DIM_M+DIM_M*(DIM_M-1)/2) 316 317 ! coeff of restitution input in one D array, solid wall 318 DOUBLE PRECISION DES_EN_WALL_INPUT(DIM_M) 319 DOUBLE PRECISION DES_ET_WALL_INPUT(DIM_M) 320 321 ! Hertzian collision model: 322 DOUBLE PRECISION :: E_YOUNG(DIM_M), Ew_YOUNG 323 DOUBLE PRECISION :: V_POISSON(DIM_M), Vw_POISSON 324 DOUBLE PRECISION :: HERT_KN(DIM_M, DIM_M), HERT_KWN(DIM_M) 325 DOUBLE PRECISION :: HERT_KT(DIM_M, DIM_M), HERT_KWT(DIM_M) 326 DOUBLE PRECISION :: G_MOD(DIM_M) 327 328 ! End particle-particle and particle-wall collision model parameters 329 !-----------------------------------------------------------------<<< 330 331 332 ! Particle attributes: radius, density, mass, moment of inertia 333 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: DES_RADIUS !(PARTICLES) 334 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: RO_Sol !(PARTICLES) 335 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: PVOL !(PARTICLES) 336 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: PMASS !(PARTICLES) 337 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: OMOI !(PARTICLES) 338 339 ! Additional quantities 340 DOUBLE PRECISION :: MIN_RADIUS, MAX_RADIUS 341 342 343 ! 'solids phase' particle diameters 344 DOUBLE PRECISION DES_D_P0 (DIM_M) 345 ! 'solids phase' particle densities 346 DOUBLE PRECISION DES_RO_s (DIM_M) 347 ! number of 'solids phases' 348 INTEGER DES_MMAX 349 350 ! Old and new particle positions, velocities (translational and 351 ! rotational) 352 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_POS_OLD !(PARTICLES,3) 353 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_POS_NEW !(PARTICLES,3) 354 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_VEL_OLD !(PARTICLES,3) 355 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_VEL_NEW !(PARTICLES,3) 356 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: OMEGA_OLD !(PARTICLES,3) 357 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: OMEGA_NEW !(PARTICLES,3) 358 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: PPOS !(PARTICLES,3) 359 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_ACC_OLD !(PARTICLES,3) 360 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: ROT_ACC_OLD !(PARTICLES,3) 361 362 ! Particle orientation 363 LOGICAL :: PARTICLE_ORIENTATION = .FALSE. 364 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: ORIENTATION !(3,PARTICLES) 365 DOUBLE PRECISION, DIMENSION(3) :: INIT_ORIENTATION = (/0.0,1.0,0.0/) 366 367 ! Defining user defined allocatable array 368 INTEGER :: DES_USR_VAR_SIZE 369 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_USR_VAR !(PARTICLES,3) 370 371 ! Total force and torque on each particle 372 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: FC !(3,PARTICLES) 373 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: TOW !(3,PARTICLES) 374 375 ! Dynamic information related to computational (eulerian) fluid grid 376 !----------------------------------------------------------------->>> 377 ! Dynamic variable. for each ijk computational fluid cell store the 378 ! total number of particles and the id's of the particles in that cell 379 TYPE iap1 380 INTEGER, DIMENSION(:), POINTER:: p 381 END TYPE iap1 382 383 ! particle can collide with at most COLLISION_ARRAY_MAX facets simultaneously 384 INTEGER, PARAMETER :: COLLISION_ARRAY_MAX = 8 385 386 ! -1 value indicates no collision 387 INTEGER, DIMENSION(:,:), ALLOCATABLE :: wall_collision_facet_id ! (COLLISION_ARRAY_MAX,PARTICLES) 388 DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: wall_collision_PFT ! (DIMN,COLLISION_ARRAY_MAX,PARTICLES) 389 390 TYPE cnaa1 391 INTEGER, DIMENSION(:), ALLOCATABLE:: p 392 INTEGER, DIMENSION(:), ALLOCATABLE:: extentdir 393 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: extentmin 394 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: extentmax 395 END TYPE cnaa1 396 ! in order to facilitate the parallel processing the PIC is defined 397 ! as single array IJK 398 TYPE(iap1), DIMENSION(:), ALLOCATABLE:: pic ! (DIMENSION_3) 399 400 TYPE(cnaa1), DIMENSION(:), ALLOCATABLE :: CELLNEIGHBOR_FACET 401 402 ! Store the number of particles in a computational fluid cell 403 INTEGER, DIMENSION(:), ALLOCATABLE :: PINC ! (DIMENSION_3) 404 405 ! For each particle track its i, j, k & ijk location on the fluid grid 406 ! and solids phase no.: 407 INTEGER, DIMENSION(:,:), ALLOCATABLE :: PIJK ! (PARTICLES,5)=>I,J,K,IJK,M 408 !-----------------------------------------------------------------<<< 409 410 411 ! note that thse variables are needed since the existing variables (i.e. 412 ! f_gs, f_ss, etc) are also being used to store the information between 413 ! the gas and continuous solids phases. 414 ! drag coefficient between gas phase and discrete particle 'phases' 415 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: F_GDS 416 ! drag coefficient between continuous solids phases and discrete 417 ! particle 'phases' 418 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: F_SDS 419 420 ! the following should probably be local to the subroutine 421 ! solve_vel_star they are only needed when invoking the non-interpolated 422 ! version of drag wherein they are used to determine part of the source 423 ! in the gas-phase momentum balances, in the routine to determine 424 ! coefficients for the pressure correction equation and in the partial 425 ! elimination algorithm 426 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: VXF_GDS 427 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: VXF_SDS 428 429 ! the contribution of solids-particle drag to the mth phase continuum 430 ! solids momentum A matrix (center coefficient) 431 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: SDRAG_AM 432 ! the contribution of solids-particle drag to the to mth phase continuum 433 ! solids momentum B vector 434 DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: SDRAG_BM 435 436 437 ! the coefficient add to gas momentum A matrix 438 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: DRAG_AM 439 ! the coefficient add to gas momentum B matrix 440 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DRAG_BM 441 442 ! Explicitly calculated fluid-particle drag force. 443 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DRAG_FC 444 445 ! An intermediate array used in calculation of mean solids velocity 446 ! by backward interpolation, i.e., when INTERP_DES_MEAN_FIELDS is true. 447 DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE ::DES_VEL_NODE 448 !(DIMENSION_3,3,DES_MMAX) 449 450 ! An intermediate array used in calculation of solids volume fraction 451 ! by backward interpolation, i.e., when INTERP_DES_MEAN_FIELDS is true. 452 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_ROPS_NODE 453 !(DIMENSION_3,DES_MMAX) 454 455 DOUBLE PRECISION, DIMENSION(:,:,:), POINTER :: weightp 456 457 ! the gas-particle drag coefficient 458 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: f_gp 459 !(PARTICLES) 460 DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: wtderivp 461 DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: sstencil 462 DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: gstencil, vstencil, pgradstencil 463 464 ! stencil for interpolation of solids pressure 465 DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: psgradstencil 466 467 ! stencil for interpolation of solids velocity 468 DOUBLE PRECISION, DIMENSION(:,:,:,:,:), ALLOCATABLE:: VEL_SOL_STENCIL 469 470 471 ! quantities are set in subroutine set_interpolation_scheme 472 ! order = order of the interpolation method, ob2l = (order+1)/2, 473 ! ob2r = order/2 474 CHARACTER(LEN=7):: scheme, interp_scheme 475 INTEGER:: order, ob2l, ob2r 476 477 ! END of interpolation related data 478 !-----------------------------------------------------------------<<< 479 480 ! Volume of each node. Used to obtain Eulerian fields 481 double precision, allocatable, dimension(:) :: des_vol_node 482 !(DIMENSION_3) 483 484 ! Ratio of actual volume of each node to volume of node not corrected for 485 ! outside the domain or being in cut-cell 486 double precision, allocatable, dimension(:) :: des_vol_node_ratio 487 !(DIMENSION_3) 488 489 ! Variable to track pressure force in computational fluid cell (ijk) 490 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: P_FORCE 491 !(DIMENSION_3,DES_MMAX) 492 493 ! Bulk density of particles in fluid cell 494 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_ROP_S,& 495 DES_ROP_SO 496 !(DIMENSION_3,DES_MMAX) 497 498 ! Volume averaged solids velocity in a fluid cell 499 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_U_s 500 !(DIMENSION_3,DES_MMAX) 501 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_V_s 502 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_W_s 503 504 ! Granular temperature in a fluid cell 505 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_THETA 506 !(DIMENSION_3,DES_MMAX) 507 508 ! Global average velocity: obtained by averaging over all the particles 509 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: DES_VEL_AVG 510 !(3) 511 512 ! Global granular energy & temp: obtained by averaging over all the particles 513 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: GLOBAL_GRAN_ENERGY 514 !(3) 515 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: GLOBAL_GRAN_TEMP 516 !(3) 517 518 ! Kinetic and potential energy of the system: obtained by averaging 519 ! over all particles 520 ! Added rotational kinetic energy (DES_ROTE) 521 DOUBLE PRECISION DES_KE, DES_PE, DES_ROTE 522 523 ! Logic for bed height calculations (T = turn on bed height 524 ! calculations) 525 LOGICAL DES_CALC_BEDHEIGHT 526 ! Used to track bed height of solids phase M 527 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: bed_height 528 !(DES_MMAX) 529 530 ! MAX velocity of particles in each direction 531 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: DES_VEL_MAX 532 !(3) 533 534 !flag to convert the outside box to facets if stl facet representation is used 535 !default is false 536 537 Logical :: des_convert_box_to_facets 538 539 Integer, parameter :: FACET_TYPE_NORMAL = 1, FACET_TYPE_PO = 2 & 540 , FACET_TYPE_MI = 3 541 542 !make this a short integer 543 !array to specify the facet type 544 Integer, dimension(:), allocatable :: STL_FACET_TYPE 545 546 Integer :: count_facet_type_normal, count_facet_type_po, count_facet_type_mi 547 548 ! Flag to turn on/off optimizing the list of facets at each des grid cell 549 550 LOGICAL :: MINIMIZE_DES_FACET_LIST 551 !-----------------------------------------------------------------<<< 552 553 554 ! Start Cohesion 555 !----------------------------------------------------------------->>> 556 ! Includes square-well type model and a van der waals type model 557 558 ! Switch to turn cohesion on and off (set in mfix.dat) 559 LOGICAL USE_COHESION 560 561 562 ! Van der Waals constants (set in mfix.dat) 563 LOGICAL SQUARE_WELL, VAN_DER_WAALS 564 DOUBLE PRECISION HAMAKER_CONSTANT 565 DOUBLE PRECISION VDW_INNER_CUTOFF ! (in cm) 566 DOUBLE PRECISION VDW_OUTER_CUTOFF 567 DOUBLE PRECISION WALL_HAMAKER_CONSTANT 568 DOUBLE PRECISION WALL_VDW_INNER_CUTOFF 569 DOUBLE PRECISION WALL_VDW_OUTER_CUTOFF 570 DOUBLE PRECISION Asperities ! average radius of asperities (default zero) 571 572 ! Store postcohesive 573 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: PostCohesive 574 !(PARTICLES) 575 ! Store cluster information array for postprocessing 576 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: PostCluster 577 578 ! Variables for van der waals cohesion calculations: 579 ! Surface energy used to calculate cohesive force for low separation distances 580 ! in Van der Waals model (this variable is calculated at the beginning of each 581 ! simulation to ensure the van der Waals force is continuous at the inner cutoff) 582 DOUBLE PRECISION SURFACE_ENERGY 583 DOUBLE PRECISION WALL_SURFACE_ENERGY 584 585 ! END Cohesion 586 !-----------------------------------------------------------------<<< 587 588 LOGICAL :: DES_EXPLICITLY_COUPLED 589 590 ! particle in cell related variable 591 type iap2 592 integer :: isize 593 integer, dimension(:), pointer:: p 594 end type iap2 595 596 type(iap2), dimension(:),allocatable:: dg_pic 597 integer, dimension(:),allocatable :: dg_pijk,dg_pijkprv 598 599 ! variable to clean the ghost cells 600 logical,dimension(:),allocatable :: ighost_updated 601 integer :: max_isize 602 603 CONTAINS 604 605 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv 606 ! ! 607 ! Subroutine: DES_CROSSPRDCT ! 608 ! Purpose: Calculate the cross product of two 3D vectors. ! 609 ! ! 610 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 611 FUNCTION DES_CROSSPRDCT(XX,YY) 612 613 IMPLICIT NONE 614 615 ! Dummy arguments 616 !---------------------------------------------------------------------// 617 ! Input vectors 618 DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: XX, YY 619 ! Result: cross product of vectors 620 DOUBLE PRECISION, DIMENSION(3) :: DES_CROSSPRDCT 621 !......................................................................! 622 623 DES_CROSSPRDCT(1) = XX(2)*YY(3) - XX(3)*YY(2) 624 DES_CROSSPRDCT(2) = XX(3)*YY(1) - XX(1)*YY(3) 625 DES_CROSSPRDCT(3) = XX(1)*YY(2) - XX(2)*YY(1) 626 627 END FUNCTION DES_CROSSPRDCT 628 629 END MODULE DISCRETELEMENT 630