1 ! -*- f90 -*- 2 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 3 ! C 4 ! Module name: DISCRETELEMENT C 5 ! Purpose: DES mod file C 6 ! Common Block containing DEM conditions C 7 ! C 8 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 9 10 MODULE DISCRETELEMENT 11 12 !----------------------------------------------- 13 ! Modules 14 !----------------------------------------------- 15 USE param, only: dim_m 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 INTEGER, PARAMETER :: DIM_M_TRI = 55 ! 10th triangular number 33 34 ! Total number of particles in simulation: read from input or generated 35 INTEGER PARTICLES 36 37 ! Start particle tracking quantities 38 !----------------------------------------------------------------->>> 39 ! Generally for inlet/outlet related routines but also employed in 40 ! tracking for parallelization 41 42 ! Dynamic particle states: 43 44 ! NONEXISTENT: This state is used with the inlet/outlet to skip 45 ! indices that do not represent particles in the system or indices 46 ! that represent particles that have exited the system. 47 48 ! ENTERING: This state identifies a particle as 'new' if true. 49 ! Particles with a classification of 'new' do not react when in contact 50 ! with a wall or another particle, however existing particles do collide 51 ! and interact with 'new' particles. The classification allows new 52 ! particles to push particles already in the system out of the way when 53 ! entering to prevent overlap. This flag is also used when the center 54 ! of a particle crosses a dem outlet (i.e. an exiting particle; see 55 ! EXITING) so that the particle will maintain its present trajectory 56 ! until it has fully exited the system 57 58 ! EXITING: This state identifies a particle as 'exiting' if true. 59 ! If a particle initiates contact with a wall surface designated as a 60 ! des outlet, this flag is set to true. With this classification the 61 ! location of the particle is checked to assess if the particle has 62 ! fully exited the system. At this point, the particle is removed 63 ! from the list. 64 65 ! GHOST, ENTERING_GHOST, EXITING_GHOST: for ghost particles 66 67 INTEGER(KIND=1), DIMENSION(:), ALLOCATABLE :: PARTICLE_STATE ! (PARTICLES) 68 69 INTEGER, PARAMETER :: nonexistent=0 70 INTEGER, PARAMETER :: normal_particle=1 71 INTEGER, PARAMETER :: entering_particle=2 72 INTEGER, PARAMETER :: exiting_particle=3 73 INTEGER, PARAMETER :: normal_ghost=4 74 INTEGER, PARAMETER :: entering_ghost=5 75 INTEGER, PARAMETER :: exiting_ghost=6 76 77 ! PARALLEL PROCESSING: explanation of variables in parallel architecture 78 ! pip - particles in each processor (includes the ghost particles) 79 ! max_pip - maximum allocated particles in processor 80 81 ! Number of particles in the system (current) 82 INTEGER PIP 83 ! Global sum of particles (excluding ghost) in the system 84 INTEGER TOT_PAR 85 ! Maximum particles permitted in the system at once 86 INTEGER MAX_PIP 87 88 ! End particle tracking quantities 89 !-----------------------------------------------------------------<<< 90 91 ! For parallel processing: global id of particles 92 INTEGER, DIMENSION(:), ALLOCATABLE :: IGLOBAL_ID 93 ! Ghost count of particles on each processor 94 INTEGER :: IGHOST_CNT 95 ! Maximum global id, new particles global id will be assigned based 96 ! on this value 97 Integer :: imax_global_id 98 99 100 ! If gener_part_config is true, then the particle_input.dat file does 101 ! not need to be supplied nor does the total number of particles as 102 ! these are determined based on the specified volume fraction (vol_frac) 103 ! in the specified domain (des_eps_xyzstart) 104 LOGICAL :: GENER_PART_CONFIG 105 DOUBLE PRECISION :: VOL_FRAC(DIM_M) 106 DOUBLE PRECISION :: DES_EPS_XSTART, & 107 DES_EPS_YSTART, DES_EPS_ZSTART 108 ! volume of the IC region for computing number of particles to be seeded 109 DOUBLE PRECISION, dimension(:), allocatable :: VOL_IC_REGION!(DIMENSION_IC) 110 ! number of particles for each phase corresponding to the IC number. This will 111 ! be real particles for DEM but parcels or computational particles for PIC model 112 INTEGER, dimension(:,:), allocatable :: PART_MPHASE_BYIC!(DIMENSION_IC, DIM_M) 113 114 ! Number of real particles by IC and by solid phase. Only relevant for PIC model 115 double precision, dimension(:,:), allocatable :: REALPART_MPHASE_BYIC!(DIMENSION_IC, DIM_M) 116 ! The number of particles that belong to solid phase M according to the 117 ! vol_frac and particle diameter. this information is used when 118 ! gener_part_config is invoked for initialization 119 ! This will be removed soon as PART_MPHASE_BYIC will be used from now on 120 INTEGER PART_MPHASE(DIM_M) 121 122 ! The number of real particles that belong to solid phase M during the initialization. 123 ! It is equal to Part_mphase for DEM but implies real number of particles for PIC model 124 double precision REALPART_MPHASE(DIM_M) 125 ! Assigns the initial particle velocity distribution based on user 126 ! specified mean and standard deviation (regardless if already set 127 ! within particle_input.dat) 128 DOUBLE PRECISION pvel_mean, PVEL_StDev 129 130 ! Output/debug controls 131 !----------------------------------------------------------------->>> 132 ! Logic that controls whether to print data dem simulations (granular or 133 ! coupled) 134 LOGICAL PRINT_DES_DATA 135 CHARACTER(LEN=255) :: VTP_DIR 136 137 ! logic that controls if des run time messages are printed on screen or not 138 LOGICAL PRINT_DES_SCREEN 139 140 ! This specifies the file type used for outputting DES data 141 ! options are : 142 ! TECPLOT - data is written in Tecplot format 143 ! undefined - data is written in ParaView format (default) 144 CHARACTER(LEN=64) :: DES_OUTPUT_TYPE 145 146 ! Used sporadically to control screen dumps (for debug purposes) 147 LOGICAL :: DEBUG_DES 148 149 ! Single particle no. index that is followed if debugging 150 INTEGER FOCUS_PARTICLE 151 152 ! Output file count for .vtp type files and for tecplot files; 153 ! for vtp output used to label .vtp file names in sequential order 154 ! and is saved so restarts begin at the correct count 155 INTEGER VTP_FINDEX, TECPLOT_FINDEX 156 ! End Output/debug controls 157 !-----------------------------------------------------------------<<< 158 159 ! DES - Continuum 160 LOGICAL DISCRETE_ELEMENT 161 LOGICAL DES_CONTINUUM_COUPLED 162 163 ! DES - Invoke hybrid model where both the DEM and continuum model 164 ! are employed to describe solids 165 LOGICAL DES_CONTINUUM_HYBRID 166 167 ! DES - 168 ! With this logic the particles see the fluid but the fluid does 169 ! not see the particles. 170 LOGICAL DES_ONEWAY_COUPLED 171 172 ! Only used when coupled and represents the number of times a pure 173 ! granular flow simulation is run before the actual coupled simulation 174 ! is started (i.e. for particle settling w/o fluid forces) 175 INTEGER NFACTOR 176 177 ! Drag 178 LOGICAL TSUJI_DRAG 179 180 ! Collision model, options are as follows 181 ! linear spring dashpot model (default/undefined) 182 ! 'hertzian' model 183 CHARACTER(LEN=64) :: DES_COLL_MODEL 184 INTEGER DES_COLL_MODEL_ENUM 185 INTEGER,PARAMETER :: HERTZIAN=0 186 INTEGER,PARAMETER :: LSD=1 187 188 ! Integration method, options are as follows 189 ! 'euler' first-order scheme (default) 190 ! 'adams_bashforth' second-order scheme (by T.Li) 191 CHARACTER(LEN=64) :: DES_INTG_METHOD 192 LOGICAL INTG_ADAMS_BASHFORTH 193 LOGICAL INTG_EULER 194 195 ! Value of solids time step based on particle properties 196 DOUBLE PRECISION DTSOLID 197 ! Run time value of simulation time used in dem simulation 198 DOUBLE PRECISION S_TIME 199 200 201 ! Neighbor search related quantities 202 !----------------------------------------------------------------->>> 203 ! Neighbor search method, options are as follows 204 ! 1= nsquare, 2=quadtree, 3=octree, 4=grid/cell based search 205 INTEGER DES_NEIGHBOR_SEARCH 206 207 ! Quantities used to determine whether neighbor search should be called 208 INTEGER NEIGHBOR_SEARCH_N 209 DOUBLE PRECISION NEIGHBOR_SEARCH_RAD_RATIO 210 LOGICAL DO_NSEARCH 211 212 ! Flag on whether to have DES_*_OLD arrays, if either Adams Bashforth or PIC is used 213 LOGICAL DO_OLD 214 215 ! Factor muliplied by sum of radii in grid based neighbor search and 216 ! nsquare search method. increases the effective radius of a particle 217 ! for detecting particle contacts 218 DOUBLE PRECISION FACTOR_RLM 219 220 ! Stores number of neighbors based on neighbor search 221 INTEGER, DIMENSION(:), ALLOCATABLE :: NEIGHBOR_INDEX 222 INTEGER, DIMENSION(:), ALLOCATABLE :: NEIGHBOR_INDEX_OLD 223 INTEGER, DIMENSION(:), ALLOCATABLE :: NEIGHBORS 224 INTEGER, DIMENSION(:), ALLOCATABLE :: NEIGHBORS_OLD 225 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: PFT_NEIGHBOR 226 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: PFT_NEIGHBOR_OLD 227 228 INTEGER :: NEIGH_NUM 229 230 ! Quantities used for reporting: max no. neighbors and max overlap 231 ! that exists during last solid time step of dem simulation 232 DOUBLE PRECISION OVERLAP_MAX 233 234 ! The number of i, j, k divisions in the grid used to perform the 235 ! cell based neighbor search 236 INTEGER :: DESGRIDSEARCH_IMAX, DESGRIDSEARCH_JMAX, & 237 DESGRIDSEARCH_KMAX 238 239 ! End neighbor search related quantities 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_TRI) 315 DOUBLE PRECISION DES_ET_INPUT(DIM_M_TRI) 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 327 ! Realistic material properties for correcting contact areas: 328 ! Actual simulations will probably use LSD with softened spring 329 ! or hertzian with unrealistically small E_Young 330 DOUBLE PRECISION :: E_YOUNG_ACTUAL(DIM_M), Ew_YOUNG_ACTUAL 331 DOUBLE PRECISION :: V_POISSON_ACTUAL(DIM_M), Vw_POISSON_ACTUAL 332 DOUBLE PRECISION :: HERT_KN_ACTUAL(DIM_M, DIM_M), HERT_KWN_ACTUAL(DIM_M) 333 334 ! Coefficients for computing binary (Hertzian) collision time. The actual 335 ! collision time is TAU_C_Base_Actual * V^(-1/5). 336 DOUBLE PRECISION :: TAU_C_Base_Actual(DIM_M, DIM_M), TAUW_C_Base_Actual(DIM_M) 337 DOUBLE PRECISION :: TAU_C_Base_Sim(DIM_M, DIM_M), TAUW_C_Base_Sim(DIM_M) 338 339 ! End particle-particle and particle-wall collision model parameters 340 !-----------------------------------------------------------------<<< 341 342 343 ! Particle attributes: radius, density, mass, moment of inertia 344 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: DES_RADIUS !(PARTICLES) 345 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: RO_Sol !(PARTICLES) 346 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: PVOL !(PARTICLES) 347 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: PMASS !(PARTICLES) 348 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: OMOI !(PARTICLES) 349 350 ! Additional quantities 351 DOUBLE PRECISION :: MIN_RADIUS, MAX_RADIUS 352 353 ! number of discrete 'solids phases' 354 INTEGER DES_MMAX 355 356 ! Old and new particle positions, velocities (translational and 357 ! rotational) 358 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_POS_OLD !(PARTICLES,3) 359 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_POS_NEW !(PARTICLES,3) 360 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_VEL_OLD !(PARTICLES,3) 361 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_VEL_NEW !(PARTICLES,3) 362 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: OMEGA_OLD !(PARTICLES,3) 363 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: OMEGA_NEW !(PARTICLES,3) 364 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: PPOS !(PARTICLES,3) 365 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_ACC_OLD !(PARTICLES,3) 366 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: ROT_ACC_OLD !(PARTICLES,3) 367 368 ! Particle orientation 369 LOGICAL :: PARTICLE_ORIENTATION = .FALSE. 370 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: ORIENTATION !(3,PARTICLES) 371 DOUBLE PRECISION, DIMENSION(3) :: INIT_ORIENTATION = (/0.0,1.0,0.0/) 372 373 ! Defining user defined allocatable array 374 INTEGER :: DES_USR_VAR_SIZE 375 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_USR_VAR !(PARTICLES,3) 376 377 ! Total force and torque on each particle 378 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: FC !(3,PARTICLES) 379 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: TOW !(3,PARTICLES) 380 381 ! particle can collide with at most COLLISION_ARRAY_MAX facets simultaneously 382 INTEGER :: COLLISION_ARRAY_MAX = 8 383 384 ! (COLLISION_ARRAY_MAX,PARTICLES) 385 ! -1 value indicates no collision 386 INTEGER, ALLOCATABLE :: wall_collision_facet_id(:,:) 387 DOUBLE PRECISION, ALLOCATABLE :: wall_collision_PFT(:,:,:) 388 389 ! Store the number of particles in a computational fluid cell 390 INTEGER, DIMENSION(:), ALLOCATABLE :: PINC ! (DIMENSION_3) 391 392 ! For each particle track its i, j, k & ijk location on the fluid grid 393 ! and solids phase no.: 394 INTEGER, DIMENSION(:,:), ALLOCATABLE :: PIJK ! (PARTICLES,5)=>I,J,K,IJK,M 395 !-----------------------------------------------------------------<<< 396 397 398 ! note that thse variables are needed since the existing variables (i.e. 399 ! f_gs, f_ss, etc) are also being used to store the information between 400 ! the gas and continuous solids phases. 401 ! drag coefficient between gas phase and discrete particle 'phases' 402 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: F_GDS 403 ! drag coefficient between continuous solids phases and discrete 404 ! particle 'phases' 405 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: F_SDS 406 407 ! the following should probably be local to the subroutine 408 ! solve_vel_star they are only needed when invoking the non-interpolated 409 ! version of drag wherein they are used to determine part of the source 410 ! in the gas-phase momentum balances, in the routine to determine 411 ! coefficients for the pressure correction equation and in the partial 412 ! elimination algorithm 413 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: VXF_GDS 414 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: VXF_SDS 415 416 ! the contribution of solids-particle drag to the mth phase continuum 417 ! solids momentum A matrix (center coefficient) 418 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: SDRAG_AM 419 ! the contribution of solids-particle drag to the to mth phase continuum 420 ! solids momentum B vector 421 DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: SDRAG_BM 422 423 424 ! the coefficient add to gas momentum A matrix 425 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: DRAG_AM 426 ! the coefficient add to gas momentum B matrix 427 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DRAG_BM 428 429 ! Explicitly calculated fluid-particle drag force. 430 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DRAG_FC 431 432 ! An intermediate array used in calculation of mean solids velocity 433 ! by backward interpolation, i.e., when INTERP_DES_MEAN_FIELDS is true. 434 DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE ::DES_VEL_NODE 435 !(DIMENSION_3,3,DIMENSION_M) 436 437 ! An intermediate array used in calculation of solids volume fraction 438 ! by backward interpolation, i.e., when INTERP_DES_MEAN_FIELDS is true. 439 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_ROPS_NODE 440 !(DIMENSION_3,DIMENSION_M) 441 442 DOUBLE PRECISION, DIMENSION(:,:,:), POINTER :: weightp 443 444 ! the gas-particle drag coefficient 445 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: f_gp 446 !(PARTICLES) 447 DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: wtderivp 448 DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: sstencil 449 DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: gstencil, vstencil, pgradstencil 450 451 ! stencil for interpolation of solids pressure 452 DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: psgradstencil 453 454 ! stencil for interpolation of solids velocity 455 DOUBLE PRECISION, DIMENSION(:,:,:,:,:), ALLOCATABLE:: VEL_SOL_STENCIL 456 457 458 ! quantities are set in subroutine set_interpolation_scheme 459 ! order = order of the interpolation method, ob2l = (order+1)/2, 460 ! ob2r = order/2 461 CHARACTER(LEN=7):: scheme, interp_scheme 462 INTEGER:: order, ob2l, ob2r 463 464 ! END of interpolation related data 465 !-----------------------------------------------------------------<<< 466 467 ! Volume of each node. Used to obtain Eulerian fields 468 double precision, allocatable, dimension(:) :: des_vol_node 469 !(DIMENSION_3,dimn,dimension_m) 470 471 ! Variable to track pressure force in computational fluid cell (ijk) 472 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: P_FORCE 473 !(DIMN, DIMENSION_3) 474 475 ! Granular temperature in a fluid cell 476 ! Global average velocity: obtained by averaging over all the particles 477 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: DES_VEL_AVG 478 !(3) 479 480 ! Global granular energy & temp: obtained by averaging over all the particles 481 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: GLOBAL_GRAN_ENERGY 482 !(3) 483 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: GLOBAL_GRAN_TEMP 484 !(3) 485 486 ! Kinetic and potential energy of the system: obtained by averaging 487 ! over all particles 488 ! Added rotational kinetic energy (DES_ROTE) 489 DOUBLE PRECISION DES_KE, DES_PE, DES_ROTE 490 491 ! Logic for bed height calculations (T = turn on bed height 492 ! calculations) 493 LOGICAL DES_CALC_BEDHEIGHT 494 ! Used to track bed height of solids phase M 495 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: bed_height 496 !(dimension_m) 497 498 ! MAX velocity of particles in each direction 499 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: DES_VEL_MAX 500 501 502 ! Flag to turn on/off optimizing the list of facets at each des grid cell 503 504 LOGICAL :: MINIMIZE_DES_FACET_LIST 505 !-----------------------------------------------------------------<<< 506 507 508 ! Start Cohesion 509 !----------------------------------------------------------------->>> 510 ! Includes square-well type model and a van der waals type model 511 512 ! Switch to turn cohesion on and off (set in mfix.dat) 513 LOGICAL USE_COHESION 514 515 516 ! Van der Waals constants (set in mfix.dat) 517 LOGICAL SQUARE_WELL, VAN_DER_WAALS 518 DOUBLE PRECISION HAMAKER_CONSTANT 519 DOUBLE PRECISION VDW_INNER_CUTOFF ! (in cm) 520 DOUBLE PRECISION VDW_OUTER_CUTOFF 521 DOUBLE PRECISION WALL_HAMAKER_CONSTANT 522 DOUBLE PRECISION WALL_VDW_INNER_CUTOFF 523 DOUBLE PRECISION WALL_VDW_OUTER_CUTOFF 524 DOUBLE PRECISION Asperities ! average radius of asperities (default zero) 525 526 ! Store postcohesive 527 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: PostCohesive 528 !(PARTICLES) 529 ! Store cluster information array for postprocessing 530 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: PostCluster 531 532 ! Variables for van der waals cohesion calculations: 533 ! Surface energy used to calculate cohesive force for low separation distances 534 ! in Van der Waals model (this variable is calculated at the beginning of each 535 ! simulation to ensure the van der Waals force is continuous at the inner cutoff) 536 DOUBLE PRECISION SURFACE_ENERGY 537 DOUBLE PRECISION WALL_SURFACE_ENERGY 538 539 ! END Cohesion 540 !-----------------------------------------------------------------<<< 541 542 LOGICAL :: DES_EXPLICITLY_COUPLED 543 544 integer, dimension(:),allocatable :: dg_pijk,dg_pijkprv 545 546 ! variable to clean the ghost cells 547 logical,dimension(:),allocatable :: ighost_updated 548 integer :: max_isize 549 550 CONTAINS 551 552 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv 553 ! ! 554 ! Subroutine: DES_CROSSPRDCT ! 555 ! Purpose: Calculate the cross product of two 3D vectors. ! 556 ! ! 557 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 558 FUNCTION DES_CROSSPRDCT(XX,YY) 559 560 IMPLICIT NONE 561 562 ! Dummy arguments 563 !---------------------------------------------------------------------// 564 ! Input vectors 565 DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: XX, YY 566 ! Result: cross product of vectors 567 DOUBLE PRECISION, DIMENSION(3) :: DES_CROSSPRDCT 568 !......................................................................! 569 570 DES_CROSSPRDCT(1) = XX(2)*YY(3) - XX(3)*YY(2) 571 DES_CROSSPRDCT(2) = XX(3)*YY(1) - XX(1)*YY(3) 572 DES_CROSSPRDCT(3) = XX(1)*YY(2) - XX(2)*YY(1) 573 574 END FUNCTION DES_CROSSPRDCT 575 576 END MODULE DISCRETELEMENT 577