calc_force_dem.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CALC_FORCE_DEM                                          !
!  Author: Jay Boyalakuntla                           Date: 12-Jun-04  !
!                                                                      !
!  Purpose: Calculate contact force and torque on particle from        !
!           particle-particle and particle-wall collisions. Treats     !
!           wall interaction also as a two-particle interaction but    !
!           accounting for the wall properties                         !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
#include "version.inc"
      SUBROUTINE calc_force_dem

! Modules
!---------------------------------------------------------------------//
      USE calc_collision_wall
      USE constant, ONLY: pi
      USE derived_types, only: multisap, boxhandle
      USE des_thermo
      USE des_thermo_cond
      USE discretelement
      USE run
      use pair_manager
      use param1, only: one, small_number, zero

      IMPLICIT NONE

! Local variables
!---------------------------------------------------------------------//
! percent of particle radius when excess overlap will be flagged
      DOUBLE PRECISION, PARAMETER :: flag_overlap = 0.20d0
! particle no. indices
      INTEGER :: I, LL, cc, CC_START, CC_END
! the overlap occuring between particle-particle or particle-wall
! collision in the normal direction
      DOUBLE PRECISION :: OVERLAP_N, OVERLAP_T(3)
! square root of the overlap
      DOUBLE PRECISION :: SQRT_OVERLAP
! distance vector between two particle centers or between a particle
! center and wall when the two surfaces are just at contact (i.e. no
! overlap)
      DOUBLE PRECISION :: R_LM,DIST_CI,DIST_CL
! the normal and tangential components of the translational relative
! velocity
      DOUBLE PRECISION :: V_REL_TRANS_NORM, rad
! distance vector between two particle centers or between a particle
! center and wall at current and previous time steps
      DOUBLE PRECISION :: DIST(3), NORMAL(3), DIST_MAG, POS(3)
! tangent to the plane of contact at current time step
      DOUBLE PRECISION :: VREL_T(3)
! normal and tangential forces
      DOUBLE PRECISION :: FN(3), FT(3)
! temporary storage of force
      DOUBLE PRECISION :: FC_TMP(3)
! temporary storage of force for torque
      DOUBLE PRECISION :: TOW_FORCE(3)
! temporary storage of torque
      DOUBLE PRECISION :: TOW_TMP(3,2)
! temporary storage of conduction/radiation
      DOUBLE PRECISION :: QQ_TMP

! store solids phase index of particle (i.e. pijk(np,5))
      INTEGER :: PHASEI, PHASELL
! local values used spring constants and damping coefficients
      DOUBLE PRECISION :: ETAN_DES, ETAT_DES
      DOUBLE PRECISION :: KN_DES, KT_DES
! local values used for calculating cohesive forces
      DOUBLE PRECISION :: FORCE_COH, EQ_RADIUS, DistApart

      LOGICAL, PARAMETER :: report_excess_overlap = .false.

      DOUBLE PRECISION :: FNMD, FTMD, MAG_OVERLAP_T, TANGENT(3)

      integer :: nn, mm, box_id, box_id2
      logical :: found
      integer :: pair(2)

!......................................................................!

! Initialize cohesive forces
      IF(use_cohesion) postcohesive(:) = zero

      CALL calc_dem_force_with_wall_stl

#ifdef do_sap
      ! do nn=0, size(multisap%saps)-1
      !    !print *,"nn = ",nn
      !    if (.not.check_boxes(multisap%saps(nn))) ERROR_STOP __LINE__
      !    if (.not.check_sort(multisap%saps(nn))) ERROR_STOP __LINE__
      ! enddo

!print *,"CALC_FORCE_DEM =================================================================================="

print *," TOTAL NUM OF NEIGHBORS IS ",neighbor_index(max_pip-1)

open (unit=123,file="neighbors.txt",action="write",status="replace")

DO ll = 1, max_pip
   cc_start = 1
   IF (ll.gt.1) cc_start = neighbor_index(ll-1)
   cc_end   = neighbor_index(ll)

   DO cc = cc_start, cc_end-1
      i  = neighbors(cc)
      write (123,*) ll,i
   enddo
enddo

close (unit=123)

             call reset_pairs(multisap%hashtable)
             do
                call get_pair(multisap%hashtable,pair)
                if (pair(1).eq.0 .and. pair(2).eq.0) exit

                if ( des_radius(pair(1))+des_radius(pair(2))> sqrt(dot_product(des_pos_new(pair(1),:)-des_pos_new(pair(2),:),des_pos_new(pair(1),:)-des_pos_new(pair(2),:)))) then
                   print *,"invalid pair: ",pair(1),pair(2)
                   stop __line__
                endif
             enddo
#endif

! Check particle LL neighbor contacts
!---------------------------------------------------------------------//

!!$omp parallel default(none) private(pos,rad,cc,cc_start,cc_end,ll,i,  &
!!$omp    overlap_n,vrel_t,v_rel_trans_norm,sqrt_overlap,dist,r_lm,     &
!!$omp    kn_des,kt_des,hert_kn,hert_kt,phasell,phasei,etan_des,        &
!!$omp    etat_des,fn,ft,overlap_t,tangent,mag_overlap_t,               &
!!$omp    eq_radius,distapart,force_coh,dist_mag,NORMAL,ftmd,fnmd,      &
!!$omp    dist_cl, dist_ci, fc_tmp, tow_tmp, tow_force, qq_tmp, box_id, box_id2, found)         &
!!$omp shared(max_pip,neighbors,neighbor_index,des_pos_new,des_radius,  &
!!$omp    des_coll_model_enum,kn,kt,pft_neighbor,pijk,                  &
!!$omp    des_etan,des_etat,mew,use_cohesion, calc_cond_des, dtsolid,   &
!!$omp    van_der_waals,vdw_outer_cutoff,vdw_inner_cutoff,              &
!!$omp    hamaker_constant,asperities,surface_energy,                   &
!!$omp    tow, fc, energy_eq, grav_mag, postcohesive, pmass, q_source, multisap, boxhandle)

!!$omp do

      DO ll = 1, max_pip
         IF(is_nonexistent(ll)) cycle
         pos = des_pos_new(ll,:)
         rad = des_radius(ll)

         cc_start = 1
         IF (ll.gt.1) cc_start = neighbor_index(ll-1)
         cc_end   = neighbor_index(ll)

         DO cc = cc_start, cc_end-1
            i  = neighbors(cc)
            IF(is_nonexistent(i)) cycle

            r_lm = rad + des_radius(i)
            dist(:) = des_pos_new(i,:) - pos(:)
            dist_mag = dot_product(dist,dist)

            fc_tmp(:) = zero

! Compute particle-particle VDW cohesive short-range forces
            IF(use_cohesion .AND. van_der_waals) THEN
               IF(dist_mag < (r_lm+vdw_outer_cutoff)**2) THEN
                  eq_radius = 2d0 * des_radius(ll)*des_radius(i) /     &
                    (des_radius(ll)+des_radius(i))
                  IF(dist_mag > (vdw_inner_cutoff+r_lm)**2) THEN
                     distapart = (sqrt(dist_mag)-r_lm)
                     force_coh = hamaker_constant * eq_radius /           &
                        (12d0*distapart**2) * (asperities/(asperities+    &
                        eq_radius) + one/(one+asperities/distapart)**2 )
                  ELSE
                     force_coh = 2d0 * pi * surface_energy * eq_radius *  &
                       (asperities/(asperities+eq_radius) + one/          &
                       (one+asperities/vdw_inner_cutoff)**2 )
                  ENDIF
                  fc_tmp(:) = dist(:)*force_coh/sqrt(dist_mag)
                  tow_tmp(:,:) = zero

! just for post-processing mag. of cohesive forces on each particle
                  postcohesive(ll) = postcohesive(ll) + force_coh / pmass(ll)
               ENDIF
            ENDIF

            IF (energy_eq) THEN
            ! Calculate conduction and radiation for thermodynamic neighbors
               IF(calc_cond_des(pijk(ll,5))) THEN
                  qq_tmp = des_conduction(ll, i, sqrt(dist_mag), pijk(ll,5), pijk(ll,4))

!!$omp atomic
                  q_source(ll) = q_source(ll) + qq_tmp

!!$omp atomic
                  q_source(i) = q_source(i) - qq_tmp
               ENDIF
            ENDIF

            IF(dist_mag > (r_lm - small_number)**2) THEN
               pft_neighbor(:,cc) = 0.0
               cycle
            ENDIF

#ifdef do_sap
               if (.not.is_pair(multisap%hashtable,ll,i)) then

                  print *,"SAP DIDNT FIND PAIR: ",ll,i
                  print *,"PARTICLE (",ll,"):  ",des_pos_new(:,ll), " WITH RADIUS: ",des_radius(ll)
                  print *,"PARTICLE (",i,"):  ",des_pos_new(:,i), " WITH RADIUS: ",des_radius(i)

                  print *,""
                  print *," ******   ",sqrt(dot_product(des_pos_new(:,ll)-des_pos_new(:,i),des_pos_new(:,ll)-des_pos_new(:,i))),"     *********"
                  print *,""

                  ! print *,"LLLLLLLLL ",boxhandle(ll)%list(:)
                  ! print *,"IIIIIIIII ",boxhandle(i)%list(:)

                  do mm=1,size(boxhandle(ll)%list)
                     if (boxhandle(ll)%list(mm)%sap_id < 0 ) cycle
                     print *," PARTICLE ",ll," IS IN ",boxhandle(ll)%list(mm)%sap_id
                     box_id = boxhandle(ll)%list(mm)%box_id

                     found = .false.
                     do nn=1,size(boxhandle(i)%list)
                        if (boxhandle(i)%list(nn)%sap_id .eq. boxhandle(ll)%list(mm)%sap_id) then
                           ! print *," PARTICLE ",i," IS ALSO IN ",boxhandle(i)%list(nn)
                           box_id2 = boxhandle(i)%list(nn)%box_id
                           found = .true.
                        endif
                     enddo

                     if (.not.found) cycle

                     ! print *,"BOTH ",ll,i," ARE IN ",boxhandle(ll)%list(mm)

                  enddo

                  error_stop __line__
               endif
#endif

            IF(dist_mag == 0) THEN
               WRITE(*,8550) ll, i
               error_stop "division by zero"
 8550 FORMAT('distance between particles is zero:',2(2x,i10))
            ENDIF

            dist_mag = sqrt(dist_mag)
            normal(:)= dist(:)/dist_mag

! Calcuate the normal overlap
            overlap_n = r_lm-dist_mag
            IF(report_excess_overlap) CALL print_excess_overlap

! Calculate the components of translational relative velocity for a
! contacting particle pair and the tangent to the plane of contact
            CALL cfrelvel(ll, i, v_rel_trans_norm, vrel_t,            &
               normal(:), dist_mag)

            phasell = pijk(ll,5)
            phasei = pijk(i,5)

! Hertz spring-dashpot contact model
            IF (des_coll_model_enum .EQ. hertzian) THEN
               sqrt_overlap = sqrt(overlap_n)
               kn_des = hert_kn(phasell,phasei)*sqrt_overlap
               kt_des = hert_kt(phasell,phasei)*sqrt_overlap
               sqrt_overlap = sqrt(sqrt_overlap)
               etan_des = des_etan(phasell,phasei)*sqrt_overlap
               etat_des = des_etat(phasell,phasei)*sqrt_overlap

! Linear spring-dashpot contact model
            ELSE
               kn_des = kn
               kt_des = kt
               etan_des = des_etan(phasell,phasei)
               etat_des = des_etat(phasell,phasei)
            ENDIF

! Calculate the normal contact force
            fn(:) =  -(kn_des * overlap_n * normal(:) + &
               etan_des * v_rel_trans_norm * normal(:))

! Calcuate the tangential overlap
            overlap_t(:) = dtsolid*vrel_t(:) + pft_neighbor(:,cc)
            mag_overlap_t = sqrt(dot_product(overlap_t,overlap_t))

! Calculate the tangential contact force.
            IF(mag_overlap_t > 0.0) THEN
! Tangential froce from spring.
               ftmd = kt_des*mag_overlap_t
! Max force before the on set of frictional slip.
               fnmd = mew*sqrt(dot_product(fn,fn))
! Direction of tangential force.
               tangent = overlap_t/mag_overlap_t
! Frictional slip
               IF(ftmd > fnmd) THEN
                  ft = -fnmd * tangent
                  overlap_t = (fnmd/kt_des) * tangent
               ELSE
                  ft = -ftmd * tangent
               ENDIF
            ELSE
               ft = 0.0
            ENDIF

! Add in the tangential dashpot damping force
            ft = ft - etat_des * vrel_t(:)

! Save tangential displacement history
            pft_neighbor(:,cc) = overlap_t(:)

! calculate the distance from the particles' centers to the contact point,
! which is taken as the radical line
! dist_ci+dist_cl=dist_li; dist_ci^2+a^2=ri^2;  dist_cl^2+a^2=rl^2
            dist_cl = dist_mag/2.d0 + (des_radius(ll)**2 - &
               des_radius(i)**2)/(2.d0*dist_mag)

            dist_ci = dist_mag - dist_cl

            tow_force(:) = des_crossprdct(normal(:), ft(:))
            tow_tmp(:,1) = dist_cl*tow_force(:)
            tow_tmp(:,2) = dist_ci*tow_force(:)

! Calculate the total force FC of a collision pair
! total contact force ( FC_TMP may already include cohesive force)
            fc_tmp(:) = fc_tmp(:) + fn(:) + ft(:)

            fc(ll,:) = fc(ll,:) + fc_tmp(:)

!!$omp atomic
            fc(i,1) = fc(i,1) - fc_tmp(1)
!!$omp atomic
            fc(i,2) = fc(i,2) - fc_tmp(2)
!!$omp atomic
            fc(i,3) = fc(i,3) - fc_tmp(3)

! for each particle the signs of norm and ft both flip, so add the same torque
            tow(ll,:) = tow(ll,:) + tow_tmp(:,1)

!!$omp atomic
            tow(i,1)  = tow(i,1)  + tow_tmp(1,2)
!!$omp atomic
            tow(i,2)  = tow(i,2)  + tow_tmp(2,2)
!!$omp atomic
            tow(i,3)  = tow(i,3)  + tow_tmp(3,2)

         ENDDO
      ENDDO
!!$omp end do

!!$omp end parallel

! just for post-processing mag. of cohesive forces on each particle
      IF(use_cohesion .AND. van_der_waals .AND. grav_mag > zero) THEN
         postcohesive(:) = postcohesive(:)/grav_mag
      ENDIF

      RETURN

      contains

        include 'functions.inc'

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: print_excess_overlap                                    !
!                                                                      !
!  Purpose: Print overlap warning messages.                            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE print_excess_overlap

      use error_manager

      IF(overlap_n > flag_overlap*des_radius(ll) .OR.                  &
         overlap_n > flag_overlap*des_radius(i)) THEN

         WRITE(err_msg,1000) trim(ival(ll)), trim(ival(i)), s_time,    &
            des_radius(ll), des_radius(i), overlap_n

         CALL flush_err_msg(header=.false., footer=.false.)
      ENDIF

 1000 FORMAT('WARNING: Excessive overplay detected between ',          &
         'particles ',a,' and ',/a,' at time ',g11.4,'.',/             &
         'RADII:  ',g11.4,' and ',g11.4,4x,'OVERLAP: ',g11.4)

      END SUBROUTINE print_excess_overlap

    END SUBROUTINE calc_force_dem