1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv 2 ! 3 ! Subroutine: CFRELVEL 4 ! Purpose: Calculate the normal and tangential components of the 5 ! relative velocity between contacting particles 6 ! 7 ! Author: Jay Boyalakuntla Date: 12-Jun-04 8 ! Reviewer: Rahul Garg Date: 01-Aug-07 9 ! 10 ! Comments: Relative (translational) velocity required for eqn 6 11 ! from the following paper: 12 ! Tsuji Y., Kawaguchi T., and Tanak T., "Lagrangian numerical 13 ! simulation of plug glow of cohesionless particles in a 14 ! horizontal pipe", Powder technology, 71, 239-250, 1992 15 ! 16 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 17 18 SUBROUTINE CFRELVEL(L, II, VRN, VSLIP, NORM, DIST_LI) 19 20 !----------------------------------------------- 21 ! Modules 22 !----------------------------------------------- 23 USE discretelement, only: DES_VEL_NEW, DES_RADIUS, OMEGA_NEW, DES_CROSSPRDCT 24 USE param1, only: ZERO 25 IMPLICIT NONE 26 !----------------------------------------------- 27 ! Dummy arguments 28 !----------------------------------------------- 29 ! indices of particle-particle contact pair 30 INTEGER, INTENT(IN) :: L, II 31 ! distance between particle centers 32 DOUBLE PRECISION, INTENT(IN) :: DIST_LI 33 ! unit normal vector along the line of contact pointing from 34 ! particle L to particle II 35 DOUBLE PRECISION, INTENT(IN) :: NORM(3) 36 ! slip velocity at point of contact 37 DOUBLE PRECISION, INTENT(OUT) :: VSLIP(3) 38 ! normal component of relative contact velocity (scalar) 39 DOUBLE PRECISION, INTENT(OUT) :: VRN 40 !----------------------------------------------- 41 ! Local variables 42 !----------------------------------------------- 43 ! translational relative velocity 44 DOUBLE PRECISION :: VRELTRANS(3) 45 ! rotational velocity at point of contact 46 DOUBLE PRECISION :: V_ROT(3), OMEGA_SUM(3) 47 ! distance from the contact point to the particle centers 48 DOUBLE PRECISION :: DIST_CL, DIST_CI 49 !----------------------------------------------- 50 51 ! translational relative velocity 52 VRELTRANS(:) = (DES_VEL_NEW(:,L) - DES_VEL_NEW(:,II)) 53 54 ! calculate the distance from the particle center to the contact point, 55 ! which is taken as the radical line 56 ! dist_ci+dist_cl=dist_li; dist_ci^2+a^2=ri^2; dist_cl^2+a^2=rl^2 57 DIST_CL = (DIST_LI**2 + DES_RADIUS(L)**2 - DES_RADIUS(II)**2)/& 58 (2.d0*DIST_LI) 59 DIST_CI = DIST_LI - DIST_CL 60 61 OMEGA_SUM(:) = OMEGA_NEW(:,L)*DIST_CL + & 62 OMEGA_NEW(:,II)*DIST_CI 63 64 ! calculate the rotational relative velocity 65 V_ROT = DES_CROSSPRDCT(OMEGA_SUM, NORM) 66 67 ! total relative velocity 68 VRELTRANS(:) = VRELTRANS(:) + V_ROT(:) 69 70 ! normal component of relative velocity (scalar) 71 VRN = DOT_PRODUCT(VRELTRANS,NORM) 72 73 ! slip velocity of the contact point 74 ! Equation (8) in Tsuji et al. 1992 75 VSLIP(:) = VRELTRANS(:) - VRN*NORM(:) 76 77 RETURN 78 END SUBROUTINE CFRELVEL 79 80 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv 81 ! 82 ! Subroutine: CFRELVEL_WALL 83 ! Purpose: Calculate the normal and tangential components of the 84 ! relative velocity between a particle and wall contact 85 ! 86 ! Comments: this subroutine is the same as above but it used for 87 ! particle-wall contact rather than particle-particle 88 ! contact. so the wall velocity is passed to this routine 89 ! rather than index of a contacting particle 90 ! 91 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 92 93 SUBROUTINE CFRELVEL_WALL(L,WALL_VEL,VRN,VSLIP,NORM,DIST_LI) 94 95 !----------------------------------------------- 96 ! Modules 97 !----------------------------------------------- 98 USE discretelement, only: DES_VEL_NEW, DES_RADIUS, OMEGA_NEW, DES_CROSSPRDCT 99 USE param1, only: ZERO 100 IMPLICIT NONE 101 !----------------------------------------------- 102 ! Dummy arguments 103 !----------------------------------------------- 104 ! indices of particle 105 INTEGER, INTENT(IN) :: L 106 ! wall velocity 107 DOUBLE PRECISION, INTENT(IN) :: WALL_VEL(3) 108 ! distance between particle centers 109 DOUBLE PRECISION, INTENT(IN) :: DIST_LI 110 ! unit normal vector along the line of contact pointing from 111 ! particle L to wall 112 DOUBLE PRECISION, INTENT(IN) :: NORM(3) 113 ! slip velocity at point of contact 114 DOUBLE PRECISION, INTENT(INOUT) :: VSLIP(3) 115 ! normal component of relative contact velocity (scalar) 116 DOUBLE PRECISION, INTENT(INOUT) :: VRN 117 !----------------------------------------------- 118 ! Local variables 119 !----------------------------------------------- 120 ! translational relative velocity 121 DOUBLE PRECISION :: VRELTRANS(3) 122 ! rotational velocity at point of contact 123 DOUBLE PRECISION :: V_ROT(3), OMEGA_SUM(3) 124 ! distance from the contact point to the particle centers 125 DOUBLE PRECISION :: DIST_CL 126 !----------------------------------------------- 127 128 ! translational relative velocity 129 VRELTRANS(:) = DES_VEL_NEW(:,L) - WALL_VEL(:) 130 131 ! calculate the distance from the particle center to the wall 132 DIST_CL = DIST_LI - DES_RADIUS(L) 133 OMEGA_SUM(:) = OMEGA_NEW(:,L)*DIST_CL 134 135 ! calculate the rotational relative velocity 136 V_ROT = DES_CROSSPRDCT(OMEGA_SUM, NORM) 137 138 ! total relative velocity 139 VRELTRANS(:) = VRELTRANS(:) + V_ROT(:) 140 141 ! normal component of relative velocity (scalar) 142 VRN = DOT_PRODUCT(VRELTRANS,NORM) 143 144 ! slip velocity of the contact point 145 ! Equation (8) in Tsuji et al. 1992 146 VSLIP(:) = VRELTRANS(:) - VRN*NORM(:) 147 148 RETURN 149 END SUBROUTINE CFRELVEL_WALL 150