drag_ss_dem_noninterp.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: DRAG_SS_DEM_NONINTERP                                   !
!                                                                      !
!  Purpose: This routine is called from the DISCRETE side to calculate !
!  the solids-solids drag force acting on each particle using cell-    !
!  center continuum solids velocity. The total contact force is also   !
!  updated to include P* for over-packing.                             !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE drag_ss_dem_noninterp

! Modules
!---------------------------------------------------------------------//
! Fluid grid loop bounds.
      use compar, only: ijkstart3, ijkend3
! Fudge factor for SS drag (c.f. Gera, 2004)
      use constant, only: segregation_slope_coefficient
! The count and a list of particles in IJK
      use discretelement, only: pinc, pijk
      use derived_types, only: pic
! Particle velocity and density
      use discretelement, only: des_vel_new, ro_sol
! Particle radius and volume.
      use discretelement, only: des_radius, pvol
! Total forces acting on particle
      use discretelement, only: fc
! Gas phase volume fraction
      use fldvar, only: ep_g
! Solids pressure
      use fldvar, only: p_star
! Diameter of continuum solids phases
      use fldvar, only: d_p
! Material (ROs) and bulk density (ROPs) of continuum solids
      use fldvar, only: rop_s, ro_s
      use functions, only: is_nonexistent, is_ghost, is_entering_ghost, is_exiting_ghost
! Function to deterine if a cell contains fluid.
      use functions, only: fluid_at
! Double precision values.
      use param1, only: zero, one
! Array sizes for solids
      use param, only: dimension_m
! Number of continuum solids phases
      use physprop, only: smax
! Flag that a continuum solids can pack
      use physprop, only: close_packed
! radial distribution function
      use rdf, only: g_0
      IMPLICIT NONE

! Local variables
!---------------------------------------------------------------------//
! Loop indices:
      INTEGER :: IJK, M, NINDX, NP, L
! average solids velocity at scalar cell center in array form
      DOUBLE PRECISION :: VELCS(3, dimension_m)
! relative velocity between solids phase m and l
      DOUBLE PRECISION :: VSLP(3), VREL
! radial distribution function between phase M and L
      DOUBLE PRECISION :: G0_ML
! Sum over all phases of ratio volume fraction over particle diameter
      DOUBLE PRECISION :: EPSoDP
! DEM particle diameter calculated from radius.
      DOUBLE PRECISION :: lDP
! solid-solid drag coefficient
      DOUBLE PRECISION :: lDss
! Intermediate calculations for volume fraction.
      DOUBLE PRECISION :: OoEPg, EPg_2
! Drag force acting on each phase.
      DOUBLE PRECISION :: D_FORCE(3)

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

! Calculate the solid-solid drag for each particle.
!---------------------------------------------------------------------//
!!$omp parallel do schedule(guided, 50) default(none)              &
!!$omp shared(IJKSTART3, IJKEND3, PINC, PIC, DES_VEL_NEW,          &
!!$omp   D_P, RO_s, ROP_s, EP_G, DES_RADIUS, RO_SOL, FC,           &
!!$omp   PVOL, SEGREGATION_SLOPE_COEFFICIENT, CLOSE_PACKED, P_STAR)&
!!$omp private(IJK, OoEPg, EPg_2, EPSoDP, NP, lDP, D_FORCE, G0_ML, &
!!$omp   VELCS, VSLP, VREL, lDss, L)
      DO ijk = ijkstart3, ijkend3

         IF(.NOT.fluid_at(ijk)) cycle
         IF(pinc(ijk) == 0) cycle

         ooepg = one/ep_g(ijk)
         epg_2 = ep_g(ijk)*ep_g(ijk)

! Calculate solids volume fraction over particle diameter.
         CALL calc_epsodp(ijk, epsodp)

! Calculate the continuum solids velocity at scalar cell center.
         CALL calc_cell_center_csolids_vel(ijk, velcs)

! Calculate the solids drag for each particle in the current cell.
         DO nindx = 1,pinc(ijk)
            np = pic(ijk)%P(nindx)
! skipping indices that do not represent particles and ghost particles
            IF(is_nonexistent(np)) cycle
            IF(is_ghost(np) .OR. is_entering_ghost(np) .OR. is_exiting_ghost(np)) cycle

! Diameter of particle (not a phase diameter).
            ldp = 2.0d0*des_radius(np)

            d_force = zero
            DO m = 1, smax

! evaluating g0 - taken from G_0.f subroutine (lebowitz form)
! this section is needed to account for all solids phases until g0 for
! multiple solids types (i.e. discrete & continuum) can be addressed
! more effectively.
!               G0_ML = OoEPg + 3.0d0*EPSoDP*D_P(IJK,M)*lDP /           &
!                  (EPg_2 *(D_P(IJK,M) + lDP))
! the calculation below allows for a more generic g0 form, however,
! the calculation above is truer to the discrete nature of the info
               l = pijk(np,5)
               g0_ml = g_0(ijk,l,m)

! Relative (slip) velocity.
               vslp = des_vel_new(np,:) - velcs(:,m)
! Relative velocity magnitude.
               vrel = sqrt(dot_product(vslp,vslp))

               CALL drag_ss_syam0(ldss, d_p(ijk,m), ldp, ro_s(ijk,m),  &
                  ro_sol(np), g0_ml, vrel)

               ldss = ldss*rop_s(ijk,m)*ro_sol(np)

! accounting for particle-particle drag due to enduring contact in a
! close-packed system
               IF(close_packed(m)) ldss = ldss +                       &
                  segregation_slope_coefficient*p_star(ijk)

! Calculating the accumulated solids-solids drag force.
               d_force(:) = d_force(:) - ldss*vslp(:)

            ENDDO ! end do loop (M=1,SMAX)

            fc(np,:) = fc(np,:) + d_force(:)*pvol(np)

         ENDDO ! END DO LOOP (NP=1,MAX_PIP)
      ENDDO ! END DO LOOP (IJK=IJKSTART3, IJKEND3)
!!$omp end parallel do

      RETURN
      END SUBROUTINE drag_ss_dem_noninterp

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: DES_DRAG_SS                                             !
!                                                                      !
!  Purpose: This subroutine is called from the CONTINUUM side. It      !
!  calculate the solids-solids drag force coefficient between          !
!  continuum and discrete solids.                                      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE drag_ss_tfm_noninterp

! Modules
!---------------------------------------------------------------------//
! Fluid grid loop bounds.
      use compar, only: ijkstart3, ijkend3
! Fudge factor for SS drag (c.f. Gera, 2004)
      use constant, only: segregation_slope_coefficient
! The count and a list of particles in IJK
      use derived_types, only: pic
      use discretelement, only: pinc, pijk, des_vol_node
! Particle velocity and density
      use discretelement, only: des_vel_new, ro_sol
! Particle radius and volume.
      use discretelement, only: des_radius
! Total forces acting on particle
      use discretelement, only: sdrag_am, sdrag_bm, f_sds
! Gas phase volume fraction
      use fldvar, only: ep_g
! Diameter of continuum solids phases
      use fldvar, only: d_p
! Material (ROs) and bulk density (ROPs) of continuum solids
      use fldvar, only: rop_s, ro_s
! Solids pressure
      use fldvar, only: p_star
! Function to deterine if a cell contains fluid.
      use functions, only: fluid_at
      use functions, only: is_nonexistent, is_ghost
      use functions, only: is_entering_ghost, is_exiting_ghost
      use geometry, only: vol
! Double precision values.
      use param1, only: zero, one
! Array sizes for solids
      use param, only: dimension_m
! Flag that a continuum solids can pack
      use physprop, only: close_packed
! Number of continuum solids phases
      use physprop, only: smax
! radial distribution function
      use rdf, only: g_0
      IMPLICIT NONE

! Local variables
!---------------------------------------------------------------------//
! Loop indices:
      INTEGER :: IJK, M, NINDX, NP, L
! average solids velocity at scalar cell center in array form
      DOUBLE PRECISION :: VELCS(3, dimension_m)
! relative velocity between solids phase m and l
      DOUBLE PRECISION :: VSLP(3), VREL
! radial distribution function between phase M and L
      DOUBLE PRECISION :: G0_ML
! Sum over all phases of ratio volume fraction over particle diameter
      DOUBLE PRECISION :: EPSoDP
! DEM particle diameter calculated from radius.
      DOUBLE PRECISION :: lDP
! solid-solid drag coefficient
      DOUBLE PRECISION :: lDss
! Intermediate calculations for volume fraction.
      DOUBLE PRECISION :: OoEPg, EPg_2
! Drag force acting on each phase.
      DOUBLE PRECISION :: lFORCE
!......................................................................!

      DO ijk = ijkstart3, ijkend3

         f_sds(ijk,:) = zero
         sdrag_am(ijk,:) = zero
         sdrag_bm(ijk,:,:) = zero

         IF(.NOT.fluid_at(ijk)) cycle
         IF(pinc(ijk) == 0) cycle

         ooepg = one/ep_g(ijk)
         epg_2 = ep_g(ijk)*ep_g(ijk)

! Calculate solids volume fraction over particle diameter.
         CALL calc_epsodp(ijk, epsodp)

! Calculate the continuum solids velocity at scalar cell center.
         CALL calc_cell_center_csolids_vel(ijk, velcs)

! Calculate the solids drag for each particle in the current cell.
         DO nindx = 1,pinc(ijk)
            np = pic(ijk)%P(nindx)
! skipping indices that do not represent particles and ghost particles
            IF(is_nonexistent(np)) cycle
            IF(is_ghost(np) .OR. is_entering_ghost(np) .OR. is_exiting_ghost(np)) cycle

! Diameter of particle (not a phase diameter).
            ldp = 2.0d0*des_radius(np)

            DO m = 1, smax

! Relative (slip) velocity.
               vslp = des_vel_new(np,:) - velcs(:,m)
! Relative velocity magnitude.
               vrel = sqrt(dot_product(vslp,vslp))

! evaluating g0 - taken from G_0.f subroutine (lebowitz form)
! this section is needed to account for all solids phases until g0 for
! multiple solids types (i.e. discrete & continuum) can be addressed
! more effectively.
!               G0_ML = OoEPg + 3.0d0*EPSoDP*D_P(IJK,M)*lDP /           &
!                  (EPg_2 *(D_P(IJK,M) + lDP))
! the calculation below allows for a more generic g0 form, however,
! the calculation above is truer to the discrete nature of the info
               l = pijk(np,5)
               g0_ml = g_0(ijk,l,m)

               CALL drag_ss_syam0(ldss, d_p(ijk,m), ldp, ro_s(ijk,m),   &
                  ro_sol(np), g0_ml, vrel)

               ldss = ldss*rop_s(ijk,m)*ro_sol(np)

! accounting for particle-particle drag due to enduring contact in a
! close-packed system
               IF(close_packed(m)) ldss = ldss +                       &
                  segregation_slope_coefficient*p_star(ijk)

! Calculating the accumulated solids-solids drag force.
               lforce = ldss/vol(l)

              sdrag_am(ijk,m) = sdrag_am(ijk,m) + lforce
              sdrag_bm(ijk,:,m) = sdrag_bm(ijk,:,m) +                  &
                 lforce*des_vel_new(np,:)

            ENDDO ! end do loop (M=1,SMAX)

         ENDDO ! END DO LOOP (NP=1,MAX_PIP)

        f_sds(ijk,:) = sdrag_am(ijk,:)

      ENDDO ! END DO LOOP (IJK=IJKSTART3, IJKEND3)

      RETURN
      END SUBROUTINE drag_ss_tfm_noninterp

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CALC_CELL_CENTER_CSOLIDS_VEL                            !
!  Author: J.Musser                                   Date: 07-NOV-14  !
!                                                                      !
!  Purpose: Calculate the scalar cell center continuum solids          !
!  velocity. This code is common to the DEM and GAS calls for non-     !
!  interpolated drag routines.                                         !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_cell_center_csolids_vel(IJK, lVELFP)

! Modules
!---------------------------------------------------------------------//
! Flags and correction factors for cut momentum cells.
      use cutcell, only: cut_u_treatment_at, theta_ue, theta_ue_bar
      use cutcell, only: cut_v_treatment_at, theta_vn, theta_vn_bar
      use cutcell, only: cut_w_treatment_at, theta_wt, theta_wt_bar
! Solids velocities.
      use fldvar, only: u_s, v_s, w_s
! Functions to average momentum to scalar cell center.
      use fun_avg, only: avg_x_e, avg_y_n, avg_z_t
! Functions to lookup adjacent cells by index.
      use functions, only: im_of, jm_of, km_of
! Flag for 3D simulatoins.
      use geometry, only: do_k
      use indices, only: i_of
! Double precision parameters.
      use param1, only: zero
! Array sizes for solids
      use param, only: dimension_m
! Number of solid phases.
      use physprop, only: smax
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Fluid cell and solids phase indices
      INTEGER, INTENT(IN) :: IJK
! Fluid velocity vector at IJK cell center.
      DOUBLE PRECISION, INTENT(OUT) :: lVELFP(3, dimension_m)

! Local variables
!---------------------------------------------------------------------//
! Phase loop counter
      INTEGER :: M
! Indices of adjacent cells
      INTEGER :: IMJK, IJMK, IJKM
!......................................................................!

! Calculate the average fluid velocity at scalar cell center.
      DO m=1,smax

         imjk = im_of(ijk)
         IF(cut_u_treatment_at(imjk)) THEN
            lvelfp(1,m) = (theta_ue_bar(imjk)*u_s(imjk,m) +            &
               theta_ue(imjk)*u_s(ijk,m))
         ELSE
            lvelfp(1,m) = avg_x_e(u_s(imjk,m),u_s(ijk,m),i_of(ijk))
         ENDIF

         ijmk = jm_of(ijk)
         IF(cut_v_treatment_at(ijmk)) THEN
            lvelfp(2,m) = (theta_vn_bar(ijmk)*v_s(ijmk,m) +            &
               theta_vn(ijmk)*v_s(ijk,m))
         ELSE
            lvelfp(2,m) = avg_y_n(v_s(ijmk,m),v_s(ijk,m))
         ENDIF

         IF(do_k) THEN
            ijkm = km_of(ijk)
            IF(cut_w_treatment_at(ijkm)) THEN
               lvelfp(3,m) = (theta_wt_bar(ijkm)*w_s(ijkm,m) +         &
                  theta_wt(ijkm)* w_s(ijk,m))
            ELSE
               lvelfp(3,m) = avg_z_t(w_s(ijkm,m),w_s(ijk,m))
            ENDIF
         ELSE
            lvelfp(3,m) = zero
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE calc_cell_center_csolids_vel

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CALC_EPSoDP                                             !
!  Author: J.Musser                                   Date: 07-NOV-14  !
!                                                                      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_epsodp(IJK, lEPSoDP)

! Global Variables:
!---------------------------------------------------------------------//
! The count and a list of particles in IJK
      use derived_types, only: pic
      use discretelement, only: pinc
! Particle volume and radius
      use discretelement, only: pvol, des_radius
! Volume of scalar cell
      use geometry, only: vol
! Continuum solids diamter
      use fldvar, only: d_p
! Function to calculate continuum solids volume fraction
      use fldvar, only: ep_s
! Number of continuum solids phases
      use physprop, only: smax

      use functions

! Global Parameters:
!---------------------------------------------------------------------//
! Double precision parameters.
      use param1, only: zero

      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Fluid cell and solids phase indices
      INTEGER, INTENT(IN) :: IJK
! Fluid velocity vector at IJK cell center.
      DOUBLE PRECISION, INTENT(OUT) :: lEPSoDP

! Local variables:
!---------------------------------------------------------------------//
! Loop counters
      INTEGER :: NP, NINDX, M

! Calculate the sum of the particle volumes divided by radius.
      lepsodp = zero
      DO nindx = 1,pinc(ijk)
         np = pic(ijk)%P(nindx)
         IF(is_nonexistent(np)) cycle
         IF(is_ghost(np) .OR. is_entering_ghost(np) .OR.&
            is_exiting_ghost(np)) cycle
         lepsodp = lepsodp + pvol(np)/des_radius(np)
      ENDDO
! Convert radius to diameter and divide by cell volume.
      lepsodp = lepsodp/(2.0d0*vol(ijk))

! Add contributions from continuum solids.
      DO m = 1, smax
         lepsodp = lepsodp + ep_s(ijk,m)/d_p(ijk,m)
      ENDDO

      RETURN
      END SUBROUTINE calc_epsodp