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